Computational Explorations of the Riemann Hypothesis: Code, Visualizations, and Research Pathways (704)

import numpy as np
from scipy.special import zeta as scipy_zeta
import matplotlib.pyplot as plt

def riemann_zeta_series(s, terms=100000):
    """
    Compute Riemann zeta function using direct series summation.
    Valid for Re(s) > 1.

    Parameters:
        s: Complex number with Re(s) > 1
        terms: Number of terms to sum

    Returns:
        Approximation of zeta(s)
    """
    result = complex(0, 0)
    for n in range(1, terms + 1):
        result += 1.0 / (n ** s)
    return result


def riemann_zeta_euler(s, terms=1000):
    """
    Compute zeta using Euler product over primes.
    Demonstrates connection between zeta and primes.

    ζ(s) = Π (1 - p^(-s))^(-1) for all primes p
    """
    def sieve_primes(limit):
        is_prime = [True] * (limit + 1)
        is_prime[0] = is_prime[1] = False
        for i in range(2, int(limit**0.5) + 1):
            if is_prime[i]:
                for j in range(i*i, limit + 1, i):
                    is_prime[j] = False
        return [i for i in range(limit + 1) if is_prime[i]]

    primes = sieve_primes(terms)
    result = complex(1, 0)

    for p in primes:
        factor = 1.0 / (1.0 - p**(-s))
        result *= factor

    return result


# Verify against known values
print("Riemann Zeta Function Computations")
print("=" * 50)
print(f"ζ(2) = {riemann_zeta_series(2):.10f}")
print(f"π²/6 = {np.pi**2/6:.10f}")
print()
print(f"ζ(4) = {riemann_zeta_series(4):.10f}")
print(f"π⁴/90 = {np.pi**4/90:.10f}")
print()
print("Euler Product Comparison:")
print(f"ζ(2) via Euler = {riemann_zeta_euler(2, 10000):.10f}")

(* Wolfram Language / Mathematica Code *)

(* Define the Riemann Zeta Function *)
RiemannZeta[s_] := Sum[1/n^s, {n, 1, Infinity}]

(* Compute specific values *)
Print["ζ(2) = ", N[Zeta[2], 20]]
Print["π²/6 = ", N[Pi^2/6, 20]]

(* Find the first 10 non-trivial zeros *)
zeros = Table[
  ZetaZero[k],
  {k, 1, 10}
];

Print["First 10 non-trivial zeros:"]
TableForm[zeros, TableHeadings -> {Range[10], {"Zero Location"}}]

(* Verify zeros lie on critical line Re(s) = 1/2 *)
realParts = Re[zeros];
Print["Real parts of zeros: ", realParts]
Print["All equal to 1/2? ", AllTrue[realParts, # == 1/2 &]]

(* Plot the zeta function along critical line *)
Plot[{
  Re[Zeta[1/2 + I*t]],
  Im[Zeta[1/2 + I*t]]
}, {t, 0, 50},
  PlotLegends -> {"Re[ζ(1/2 + it)]", "Im[ζ(1/2 + it)]"},
  PlotLabel -> "Riemann Zeta on Critical Line",
  AxesLabel -> {"t", "Value"}
]

                    COMPLEX PLANE STRUCTURE
    ═══════════════════════════════════════════════════════

    Im(s)
      ↑
      │
   40 ┤                    •  ← Zero at 1/2 + 37.586i
      │                    │
   30 ┤                    •  ← Zero at 1/2 + 32.935i
      │                    │
   20 ┤                    •  ← Zero at 1/2 + 21.022i
      │          ┌─────────┼─────────┐
   14 ┤          │ CRITICAL│ STRIP   │  ← Zero at 1/2 + 14.135i
      │          │    ↓    •    ↓    │
      │          │         │         │
    0 ┼──────────┼─────────┼─────────┼──────────→ Re(s)
      │          │         │         │
      0         0.5        1        1.5
                 ↑
           CRITICAL LINE
           (Re(s) = 1/2)

    Legend:
    • = Non-trivial zeros (all on critical line if RH is true)
    █ = Critical strip (0 < Re(s) < 1)
    │ = Critical line (Re(s) = 1/2)

import numpy as np
import matplotlib.pyplot as plt
from mpmath import zetazero, zeta, mp

# Set high precision
mp.dps = 50

def plot_zeta_zeros(num_zeros=30):
    """
    Plot the first num_zeros non-trivial zeros of zeta.
    """
    # Get zeros using mpmath
    zeros = [complex(zetazero(n)) for n in range(1, num_zeros + 1)]

    real_parts = [z.real for z in zeros]
    imag_parts = [z.imag for z in zeros]

    # Create figure
    fig, axes = plt.subplots(1, 2, figsize=(14, 6))

    # Plot 1: Zeros in complex plane
    ax1 = axes[0]
    ax1.axvline(x=0.5, color='red', linestyle='--',
                linewidth=2, label='Critical Line')
    ax1.axvspan(0, 1, alpha=0.1, color='blue', label='Critical Strip')
    ax1.scatter(real_parts, imag_parts, c='green', s=50,
                zorder=5, label='Zeta Zeros')
    ax1.set_xlabel('Re(s)')
    ax1.set_ylabel('Im(s)')
    ax1.set_title('Non-trivial Zeros of ζ(s)')
    ax1.legend()
    ax1.grid(True, alpha=0.3)
    ax1.set_xlim(-0.5, 1.5)

    # Plot 2: Distribution of imaginary parts
    ax2 = axes[1]
    ax2.barh(range(len(imag_parts)), imag_parts, color='steelblue')
    ax2.set_xlabel('Imaginary Part')
    ax2.set_ylabel('Zero Index')
    ax2.set_title('Heights of Zeta Zeros')
    ax2.grid(True, alpha=0.3)

    plt.tight_layout()
    plt.savefig('zeta_zeros.png', dpi=150)
    plt.show()

    # Print statistics
    print("\nZero Statistics:")
    print(f"  Number of zeros: {num_zeros}")
    print(f"  All Re(s) = 0.5: {all(abs(r - 0.5) < 1e-10 for r in real_parts)}")
    print(f"  Largest height: {max(imag_parts):.6f}")
    print(f"  Average spacing: {np.mean(np.diff(sorted(imag_parts))):.6f}")

# Run visualization
plot_zeta_zeros(30)

import numpy as np
from mpmath import li, log, sqrt, pi, zetazero

def prime_sieve(limit):
    """
    Sieve of Eratosthenes - find all primes up to limit.
    """
    is_prime = np.ones(limit + 1, dtype=bool)
    is_prime[0:2] = False

    for i in range(2, int(np.sqrt(limit)) + 1):
        if is_prime[i]:
            is_prime[i*i::i] = False

    return np.where(is_prime)[0]


def prime_counting_function(x):
    """
    π(x) = number of primes ≤ x
    """
    if x < 2:
        return 0
    primes = prime_sieve(int(x))
    return len(primes)


def riemann_prime_formula(x, num_zeros=100):
    """
    Riemann's explicit formula for prime counting.
    Uses the first num_zeros zeta zeros.

    R(x) = Σ (μ(n)/n) * li(x^(1/n))

    With correction from zeros:
    π(x) ≈ R(x) - Σ R(x^ρ) where ρ are zeta zeros
    """
    from mpmath import mobius, mpf

    # Main term: Riemann R function
    def R(y):
        if y <= 1:
            return mpf(0)
        result = mpf(0)
        for n in range(1, 100):
            mu = mobius(n)
            if mu != 0:
                term = (mu / n) * li(y ** (1/n))
                result += term
                if abs(term) < 1e-10:
                    break
        return result

    main_term = float(R(x))

    # Correction from zeros
    correction = 0
    for k in range(1, min(num_zeros + 1, 50)):
        rho = complex(zetazero(k))
        # Add contribution from zero and its conjugate
        term = R(x ** rho)
        correction += 2 * float(term.real)

    return main_term - correction


# Compare exact count with Riemann's formula
print("Prime Counting: Exact vs Riemann Formula")
print("=" * 55)
print(f"{'x':>10} {'π(x)':>10} {'R(x)':>12} {'Error':>10}")
print("-" * 55)

for x in [100, 1000, 10000, 100000]:
    exact = prime_counting_function(x)
    approx = riemann_prime_formula(x, 30)
    error = abs(exact - approx)
    print(f"{x:>10} {exact:>10} {approx:>12.2f} {error:>10.2f}")

(* Wolfram Language: Prime-Zeta Connection *)

(* Prime counting function *)
PrimeCountExact[x_] := PrimePi[x]

(* Riemann's R function *)
RiemannR[x_] := Sum[
  MoebiusMu[n]/n * LogIntegral[x^(1/n)],
  {n, 1, 100}
]

(* Compare approximations *)
data = Table[{
    x,
    PrimePi[x],
    N[RiemannR[x], 10],
    N[LogIntegral[x], 10],
    N[x/Log[x], 10]
  },
  {x, {100, 1000, 10000, 100000, 1000000}}
];

Print["Comparison of Prime Counting Approximations:"]
TableForm[data,
  TableHeadings -> {
    None,
    {"x", "π(x)", "R(x)", "Li(x)", "x/ln(x)"}
  }
]

(* Visualize the error term *)
Plot[{
  PrimePi[x] - LogIntegral[x],
  PrimePi[x] - RiemannR[x]
}, {x, 2, 10000},
  PlotLegends -> {"π(x) - Li(x)", "π(x) - R(x)"},
  PlotLabel -> "Error in Prime Counting Approximations",
  AxesLabel -> {"x", "Error"}
]

(* Connection to RH: Error bound *)
(* If RH is true: |π(x) - Li(x)| < (1/8π) * √x * ln(x) *)
RHBound[x_] := (1/(8*Pi)) * Sqrt[x] * Log[x]

Plot[{
  Abs[PrimePi[x] - LogIntegral[x]],
  RHBound[x]
}, {x, 100, 100000},
  PlotLegends -> {"Actual Error", "RH Bound"},
  PlotLabel -> "Prime Counting Error vs RH Bound",
  PlotScale -> "Log"
]

from mpmath import mp, zeta, zetazero, diff, fabs
import numpy as np

# High precision
mp.dps = 30

def verify_zero(k, tolerance=1e-20):
    """
    Verify that the k-th zeta zero lies on critical line.

    Returns: (zero, is_on_critical_line, |zeta(zero)|)
    """
    zero = zetazero(k)
    real_part = float(zero.real)
    zeta_value = abs(complex(zeta(zero)))

    is_on_line = abs(real_part - 0.5) < tolerance
    is_zero = zeta_value < tolerance

    return {
        'index': k,
        'zero': complex(zero),
        'real_part': real_part,
        'imaginary_part': float(zero.imag),
        'on_critical_line': is_on_line,
        'zeta_value': zeta_value,
        'is_valid_zero': is_zero
    }


def gram_points(n):
    """
    Compute Gram points - where Im(ζ(1/2 + it)) changes sign.
    Useful for locating zeros.
    """
    from mpmath import grampoint
    return [float(grampoint(k)) for k in range(n)]


# Verify first 20 zeros
print("Verification of Riemann Zeta Zeros")
print("=" * 70)
print(f"{'k':>4} {'Im(ρ)':>15} {'Re(ρ)':>10} {'|ζ(ρ)|':>15} {'Valid':>8}")
print("-" * 70)

all_valid = True
for k in range(1, 21):
    result = verify_zero(k)
    valid = "✓" if result['is_valid_zero'] and result['on_critical_line'] else "✗"
    if valid == "✗":
        all_valid = False
    print(f"{k:>4} {result['imaginary_part']:>15.6f} "
          f"{result['real_part']:>10.6f} {result['zeta_value']:>15.2e} {valid:>8}")

print("-" * 70)
print(f"All zeros verified on critical line: {'YES' if all_valid else 'NO'}")

# Gram points analysis
print("\nGram Points (zero-crossing indicators):")
grams = gram_points(10)
for i, g in enumerate(grams):
    print(f"  g_{i} = {g:.6f}")

            ZERO VERIFICATION ALGORITHM
    ════════════════════════════════════════════

    ┌─────────────────────────────────────────┐
    │  INPUT: Index k of desired zero        │
    └─────────────────────┬───────────────────┘
                          │
                          ▼
    ┌─────────────────────────────────────────┐
    │  Compute Gram points g_{k-1}, g_k      │
    │  These bracket the k-th zero           │
    └─────────────────────┬───────────────────┘
                          │
                          ▼
    ┌─────────────────────────────────────────┐
    │  Use Newton-Raphson iteration:         │
    │                                         │
    │  t_{n+1} = t_n - ζ(1/2+it_n)           │
    │                  ─────────────          │
    │                  ζ'(1/2+it_n)          │
    └─────────────────────┬───────────────────┘
                          │
                          ▼
    ┌─────────────────────────────────────────┐
    │  Verify:                               │
    │  • |ζ(ρ)| < 10^{-20}                   │
    │  • |Re(ρ) - 0.5| < 10^{-20}            │
    └─────────────────────┬───────────────────┘
                          │
                          ▼
    ┌─────────────────────────────────────────┐
    │  OUTPUT: Verified zero ρ = 1/2 + it    │
    │  with certificate of accuracy          │
    └─────────────────────────────────────────┘

(* Random Matrix Theory Connection *)
(* Montgomery's Pair Correlation Conjecture *)

(* Get normalized zero spacings *)
zeros = Table[Im[ZetaZero[k]], {k, 1, 1000}];
spacings = Differences[zeros];
meanSpacing = Mean[spacings];
normalizedSpacings = spacings / meanSpacing;

(* Compare to GUE prediction *)
GUEDensity[s_] := 1 - (Sin[Pi*s]/(Pi*s))^2

(* Plot histogram vs GUE *)
Show[
  Histogram[normalizedSpacings,
    {0, 3, 0.1},
    "PDF",
    PlotLabel -> "Zero Spacing Distribution"
  ],
  Plot[GUEDensity[s], {s, 0, 3},
    PlotStyle -> {Red, Thick},
    PlotLegends -> {"GUE Prediction"}
  ]
]

(* Compute pair correlation *)
PairCorrelation[zeros_, delta_] := Module[
  {n = Length[zeros], count = 0, avgDensity},
  avgDensity = n / (Max[zeros] - Min[zeros]);

  Do[
    If[i != j && Abs[zeros[[i]] - zeros[[j]]] < delta,
      count++
    ],
    {i, n}, {j, n}
  ];

  count / (n * 2 * delta * avgDensity)
]

(* Test Montgomery's conjecture *)
Print["Pair Correlation Analysis:"]
Table[
  {delta, PairCorrelation[zeros[[1;;200]], delta]},
  {delta, {0.1, 0.2, 0.5, 1.0}}
] // TableForm

import numpy as np
from mpmath import zetazero
from scipy import stats

def analyze_zero_spacings(num_zeros=500):
    """
    Analyze statistical properties of zeta zero spacings.
    Compare with Random Matrix Theory predictions.
    """
    # Compute zeros
    print(f"Computing {num_zeros} zeta zeros...")
    zeros = [float(zetazero(k).imag) for k in range(1, num_zeros + 1)]

    # Compute spacings
    spacings = np.diff(zeros)
    mean_spacing = np.mean(spacings)
    normalized = spacings / mean_spacing

    # Statistical measures
    print("\nSpacing Statistics:")
    print(f"  Mean spacing: {mean_spacing:.6f}")
    print(f"  Std deviation: {np.std(normalized):.6f}")
    print(f"  Skewness: {stats.skew(normalized):.6f}")
    print(f"  Kurtosis: {stats.kurtosis(normalized):.6f}")

    # Compare with theoretical predictions
    # Poisson (random): mean = 1, var = 1
    # GUE (RMT): mean = 1, var ≈ 0.178

    poisson_var = 1.0
    gue_var = 0.178  # π²/18 - 1/3

    actual_var = np.var(normalized)

    print("\nComparison with Theoretical Models:")
    print(f"  Actual variance: {actual_var:.4f}")
    print(f"  Poisson variance: {poisson_var:.4f}")
    print(f"  GUE variance: {gue_var:.4f}")

    if abs(actual_var - gue_var) < abs(actual_var - poisson_var):
        print("  → Spacings match GUE (supports RH)")
    else:
        print("  → Spacings match Poisson (unexpected)")

    return normalized


# Run analysis
spacings = analyze_zero_spacings(500)

#!/usr/bin/env python3
"""
Complete Riemann Hypothesis Analysis Suite
Combines all computational methods discussed.
"""

import numpy as np
from mpmath import mp, zeta, zetazero, li, log, pi
from scipy.special import zeta as scipy_zeta
import time

mp.dps = 30  # 30 decimal places precision


class RiemannAnalyzer:
    """Comprehensive toolkit for RH analysis."""

    def __init__(self, num_zeros=100):
        self.num_zeros = num_zeros
        self.zeros = None
        self.primes = None

    def compute_zeros(self):
        """Compute and cache zeta zeros."""
        print(f"Computing {self.num_zeros} zeros...")
        start = time.time()
        self.zeros = [complex(zetazero(k)) for k in range(1, self.num_zeros + 1)]
        elapsed = time.time() - start
        print(f"Completed in {elapsed:.2f} seconds")
        return self.zeros

    def verify_critical_line(self, tolerance=1e-10):
        """Verify all zeros lie on Re(s) = 1/2."""
        if self.zeros is None:
            self.compute_zeros()

        violations = []
        for i, z in enumerate(self.zeros):
            if abs(z.real - 0.5) > tolerance:
                violations.append((i+1, z))

        if violations:
            print(f"WARNING: {len(violations)} zeros off critical line!")
            for idx, z in violations:
                print(f"  Zero {idx}: {z}")
        else:
            print(f"✓ All {len(self.zeros)} zeros on critical line")

        return len(violations) == 0

    def compute_prime_error(self, x_max=100000):
        """Compute error in prime counting formula."""
        from sympy import primepi

        results = []
        for x in [100, 1000, 10000, x_max]:
            exact = int(primepi(x))
            li_approx = float(li(x))
            error = exact - li_approx

            # RH bound: |π(x) - Li(x)| < (1/8π)√x ln(x)
            rh_bound = (1/(8*float(pi))) * np.sqrt(x) * np.log(x)

            results.append({
                'x': x,
                'exact': exact,
                'li_approx': li_approx,
                'error': error,
                'rh_bound': rh_bound,
                'within_rh': abs(error) < rh_bound
            })

        return results

    def analyze_spacings(self):
        """Statistical analysis of zero spacings."""
        if self.zeros is None:
            self.compute_zeros()

        imag_parts = [z.imag for z in self.zeros]
        spacings = np.diff(imag_parts)
        normalized = spacings / np.mean(spacings)

        return {
            'mean': np.mean(normalized),
            'variance': np.var(normalized),
            'gue_variance': 0.178,
            'matches_gue': abs(np.var(normalized) - 0.178) < 0.05
        }


# Run complete analysis
if __name__ == "__main__":
    print("=" * 60)
    print("RIEMANN HYPOTHESIS ANALYSIS SUITE")
    print("=" * 60)

    analyzer = RiemannAnalyzer(num_zeros=100)

    print("\n[1] Computing Zeta Zeros")
    print("-" * 40)
    analyzer.compute_zeros()

    print("\n[2] Verifying Critical Line")
    print("-" * 40)
    analyzer.verify_critical_line()

    print("\n[3] Prime Counting Error Analysis")
    print("-" * 40)
    results = analyzer.compute_prime_error()
    for r in results:
        status = "✓" if r['within_rh'] else "✗"
        print(f"  x={r['x']:>7}: error={r['error']:>8.2f}, "
              f"RH bound={r['rh_bound']:>10.2f} {status}")

    print("\n[4] Zero Spacing Statistics")
    print("-" * 40)
    stats = analyzer.analyze_spacings()
    print(f"  Variance: {stats['variance']:.4f}")
    print(f"  GUE prediction: {stats['gue_variance']:.4f}")
    print(f"  Matches GUE: {'Yes ✓' if stats['matches_gue'] else 'No ✗'}")

    print("\n" + "=" * 60)
    print("ANALYSIS COMPLETE")
    print("=" * 60)