E-002-first-zeros.py

{
  "experiment": "E-002",
  "task": "0.6a",
  "description": "First 10000 nontrivial zero ordinates via mpmath.zetazero (multiprocess)",
  "script": "research/numerics/scripts/E-002-first-zeros.py",
  "parameters": {
    "n_zeros": 10000,
    "dps": 25,
    "workers": 15
  },
  "checks": {
    "strictly_monotone": true,
    "gamma_1_matches_known": true,
    "riemann_von_mangoldt_abs_error": "0.3207309991"
  },
  "gap_statistics": {
    "min_gap": "0.03769849773",
    "max_gap": "6.887314497",
    "mean_gap": "0.9864634392",
    "local_mean_spacing_at_end": "0.8536743703"
  },
  "elapsed_seconds": 0.0,
  "timestamp_utc": "2026-08-17T20:09:20.100449+00:00",
  "row_count": 10000
}

"""E-002 — Task 0.6a: first 10^4 nontrivial zero ordinates of zeta.

Computes γ_n = Im(ρ_n) for n = 1..N via mpmath.zetazero at dps 25.
Parallel: multiprocessing over contiguous n-blocks (zetazero(n) is
independent per n); each worker writes its own part file with per-row
flush (crash-safe), main merges. Resumable: existing complete part files
are skipped.

Sanity checks: strict monotonicity; γ_1 = 14.134725141734693790…;
Riemann–von Mangoldt N(T) agreement at T = γ_N.

Outputs: research/numerics/data/E-002-zeros-first-10000.csv (+ .json)
Parts:   research/numerics/data/E-002-parts/part-NNN.csv
Cached:  does not recompute if final CSV is complete (use --force).
"""
import json
import os
import sys
import time
from datetime import datetime, timezone
from multiprocessing import Pool
from pathlib import Path

import mpmath as mp
import numpy as np

REPO = Path(__file__).resolve().parents[3]
DATA = REPO / "research" / "numerics" / "data"
PARTS = DATA / "E-002-parts"
N_ZEROS = 10_000
DPS = 25
N_WORKERS = min(15, os.cpu_count() or 4)  # leave one core for the OS

CSV = DATA / "E-002-zeros-first-10000.csv"
JSON = DATA / "E-002-zeros-first-10000.json"


def block_bounds(n_blocks: int):
    bounds = []
    per = N_ZEROS // n_blocks
    lo = 1
    for b in range(n_blocks):
        hi = N_ZEROS if b == n_blocks - 1 else lo + per - 1
        bounds.append((lo, hi))
        lo = hi + 1
    return bounds


def part_path(lo: int) -> Path:
    return PARTS / f"part-{lo:05d}.csv"


def compute_block(lo: int, hi: int) -> str:
    mp.mp.dps = DPS
    p = part_path(lo)
    tmp = p.with_suffix(".tmp")
    with tmp.open("w", encoding="utf-8", newline="") as f:
        for n in range(lo, hi + 1):
            g = mp.im(mp.zetazero(n))
            f.write(f"{n},{mp.nstr(g, DPS)}\n")
            f.flush()  # crash-safe: no buffered loss on kill
    os.replace(tmp, p)  # atomic: complete parts only
    return str(p)


def part_complete(lo: int, hi: int) -> bool:
    p = part_path(lo)
    if not p.exists():
        return False
    with p.open(encoding="utf-8") as f:
        return sum(1 for _ in f) == hi - lo + 1


def main() -> int:
    force = "--force" in sys.argv
    if CSV.exists() and JSON.exists() and not force:
        with CSV.open(encoding="utf-8") as f:
            if sum(1 for _ in f) - 1 >= N_ZEROS:
                print(f"cached: {CSV} complete (use --force to recompute)")
                return 0
    DATA.mkdir(parents=True, exist_ok=True)
    PARTS.mkdir(parents=True, exist_ok=True)

    blocks = block_bounds(N_WORKERS)
    todo = [(lo, hi) for lo, hi in blocks
            if force or not part_complete(lo, hi)]
    print(f"blocks: {len(blocks)} total, {len(todo)} to compute "
          f"({N_WORKERS} workers)")

    t0 = time.time()
    if todo:
        with Pool(N_WORKERS) as pool:
            for i, path in enumerate(pool.starmap(compute_block, todo), 1):
                print(f"  part done {i}/{len(todo)}: {Path(path).name} "
                      f"({time.time()-t0:.0f}s)", flush=True)

    # merge
    rows = []
    for lo, hi in blocks:
        with part_path(lo).open(encoding="utf-8") as f:
            rows.extend(f.read().strip().splitlines())
    rows.sort(key=lambda r: int(r.split(",", 1)[0]))
    with CSV.open("w", encoding="utf-8", newline="") as f:
        f.write("n,gamma_n\n")
        f.write("\n".join(rows) + "\n")
    elapsed = time.time() - t0

    # sanity checks
    mp.mp.dps = DPS
    # full-precision values straight from the CSV text (not float64 reload)
    with CSV.open(encoding="utf-8") as f:
        next(f)
        gammas = [mp.mpf(line.split(",", 1)[1]) for line in f]
    mono = all(gammas[i] < gammas[i + 1] for i in range(len(gammas) - 1))
    g1_known = mp.mpf("14.134725141734693790457251983562470270784257115699")
    g1_ok = abs(gammas[0] - g1_known) < mp.mpf("1e-20")
    T = gammas[-1]
    rvm = (T / (2 * mp.pi)) * (mp.log(T / (2 * mp.pi)) - 1) + mp.mpf("0.875")
    rvm_err = abs(N_ZEROS - rvm)
    gaps = [gammas[i + 1] - gammas[i] for i in range(len(gammas) - 1)]
    stats = {
        "min_gap": mp.nstr(min(gaps), 10),
        "max_gap": mp.nstr(max(gaps), 10),
        "mean_gap": mp.nstr(sum(gaps) / len(gaps), 10),
        "local_mean_spacing_at_end": mp.nstr(2 * mp.pi / mp.log(T / (2 * mp.pi)), 10),
    }
    print(f"monotone: {mono} | g1 ok: {g1_ok} | RvM err: "
          f"{mp.nstr(rvm_err, 5)} (should be << log T = "
          f"{mp.nstr(mp.log(T), 3)})")
    print("gap stats:", stats)

    meta = {
        "experiment": "E-002",
        "task": "0.6a",
        "description": f"First {N_ZEROS} nontrivial zero ordinates via "
                       "mpmath.zetazero (multiprocess)",
        "script": "research/numerics/scripts/E-002-first-zeros.py",
        "parameters": {"n_zeros": N_ZEROS, "dps": DPS,
                       "workers": N_WORKERS},
        "checks": {
            "strictly_monotone": bool(mono),
            "gamma_1_matches_known": bool(g1_ok),
            "riemann_von_mangoldt_abs_error": mp.nstr(rvm_err, 10),
        },
        "gap_statistics": stats,
        "elapsed_seconds": round(elapsed, 1),
        "timestamp_utc": datetime.now(timezone.utc).isoformat(),
        "row_count": N_ZEROS,
    }
    JSON.write_text(json.dumps(meta, indent=2), encoding="utf-8")
    ok = mono and g1_ok and rvm_err < 2 * mp.log(T)
    print(f"wrote {CSV.name}, {JSON.name} in {elapsed:.0f}s — "
          f"{'ALL CHECKS PASS' if ok else 'CHECK FAILURE'}")
    return 0 if ok else 1


if __name__ == "__main__":
    raise SystemExit(main())