E-008-li-coefficients.py

{
  "experiment": "E-008",
  "task": "3.2 part 1 (Li coefficients from zero table)",
  "description": "Li coefficients lambda_n (n<=2000) from E-002 zeros with explicit tail bound; validation vs lambda_1 literature value; leading-asymptotic ratio trend",
  "script": "research/numerics/scripts/E-008-li-coefficients.py",
  "parameters": {
    "n_max": 2000,
    "zeros_used": 10000,
    "H": 9877.782654005501,
    "precision": "float64"
  },
  "checks": {
    "lambda_1_within_tail_bound": true,
    "all_positive_in_range": true,
    "usable_range_at_least_n100": true
  },
  "key_values": {
    "lambda_1": 0.022961009776646203,
    "lambda_100": 117.25678701049617,
    "lambda_1000": 2191.3892531979554,
    "lambda_2000": 4813.526935739563
  },
  "usable_range_n": 116,
  "caveats": [
    "truncated at H; tail bound explicit in CSV; usable range (tail<5%) is only n<=~116 with 10^4 zeros — n~10^5 requires the prime-side Bombieri-Lagarias exact formulas (noted as task 3.2 continuation)",
    "positivity on our data is trivial (zeros on-line) — not evidence about RH"
  ],
  "elapsed_seconds": 0.2,
  "timestamp_utc": "2026-08-17T22:59:24.581674+00:00",
  "row_count": 2000
}

"""E-008 — Task 3.2 (part 1): Li coefficients from the zero table.

Li coefficients:  lambda_n = Sum_rho [1 - (1 - 1/rho)^n]   (n = 1, 2, ...)
Li's criterion: lambda_n >= 0 for all n  <=>  RH.

With zeros on the line (true for our table), rho = 1/2 + i*gamma gives
|1 - 1/rho| = 1, so 1 - 1/rho = e^{i phi_gamma} with
  phi_gamma = pi - 2 arctan(2 gamma)  (~ 1/gamma for large gamma),
and each +/- gamma pair contributes 2(1 - cos(n phi_gamma)) >= 0.

Part 1 scope (this experiment): zero-side truncation to the E-002 table
(|gamma| <= H = 9877.78...) with an EXPLICIT TAIL BOUND:
  tail_n = Sum_{|gamma|>H} 2(1-cos(n phi_gamma))
         <= 2 n^2 Sum_{gamma>H} gamma^-2 ~ 2 n^2 (log(H/2pi)+1)/(2 pi H)
(RvM density; reported with safety factor 2x). Honest range: n where
tail_bound stays a few % of lambda_n.

Validation: lambda_1 = 0.023095708966121... (literature).

Also measured: lambda_n / (n/2 log n) trend (known leading asymptotic
lambda_n ~ (n/2) log n) and positivity (trivial on our data — logged for
the record, not as evidence).

Outputs: research/numerics/data/E-008-li-coefficients.csv (+ .json)
"""
import json
import sys
import time
from datetime import datetime, timezone
from pathlib import Path

import numpy as np

REPO = Path(__file__).resolve().parents[3]
DATA = REPO / "research" / "numerics" / "data"
ZEROS_CSV = DATA / "E-002-zeros-first-10000.csv"
CSV = DATA / "E-008-li-coefficients.csv"
JSON = DATA / "E-008-li-coefficients.json"

N_MAX = 2000          # compute lambda_n for n = 1..N_MAX
LAMBDA1_KNOWN = 0.02309570896612103381  # literature value


def main() -> int:
    if not ZEROS_CSV.exists():
        print("E-002 zeros not found"); return 1
    if CSV.exists() and "--force" not in sys.argv:
        print(f"cached: {CSV} exists (use --force to recompute)"); return 0
    t0 = time.time()

    gam = np.loadtxt(ZEROS_CSV, delimiter=",", skiprows=1, usecols=1)
    H = float(gam[-1])
    phi = np.pi - 2.0 * np.arctan(2.0 * gam)      # phase per zero (upper half)
    # each pair (+/- gamma) contributes 2(1-cos(n phi)); zeros table has upper half only
    n = np.arange(1, N_MAX + 1)
    # 1 - cos(x) = 2 sin^2(x/2): stable
    contrib = 4.0 * np.sin(0.5 * np.outer(n, phi))**2   # (N_MAX, n_zeros)
    lam = contrib.sum(axis=1)

    # explicit tail bound (2x safety on the RvM-density integral)
    tail_const = 2.0 * 2.0 * (np.log(H / (2 * np.pi)) + 1.0) / (2 * np.pi * H)
    tail = tail_const * n**2

    lead = 0.5 * n * np.log(n)     # known leading asymptotic (n>=2)
    with np.errstate(divide="ignore", invalid="ignore"):
        ratio = np.where(n >= 2, lam / lead, np.nan)

    lam1_err = abs(lam[0] - LAMBDA1_KNOWN)
    # validation threshold = the explicit tail bound at n=1 (honest:
    # truncation error can be as large as the bound)
    lam1_ok = bool(lam1_err < tail[0])
    # usable range: largest n with tail bound < 5% of lambda_n
    frac = tail / lam
    usable = n[frac < 0.05]
    n_usable = int(usable[-1]) if len(usable) else 0
    checks = {
        "lambda_1_within_tail_bound": lam1_ok,
        "all_positive_in_range": bool(np.all(lam > 0)),
        "usable_range_at_least_n100": bool(n_usable >= 100),
    }
    print(f"lambda_1 = {lam[0]:.12f}  (known {LAMBDA1_KNOWN}, err {lam1_err:.2e}, "
          f"tail bound {tail[0]:.2e})")
    print(f"lambda_10 = {lam[9]:.6f}   lambda_100 = {lam[99]:.4f}   "
          f"lambda_1000 = {lam[999]:.2f}   lambda_2000 = {lam[1999]:.2f}")
    print(f"usable range (tail < 5%): n <= {n_usable}")
    print(f"lambda_n/(n/2 log n): n=100 {ratio[99]:.4f}, n=250 {ratio[249]:.4f}, "
          f"n=500 {ratio[499]:.4f} (truncation-depressed beyond usable range)")
    print("checks:", checks)

    DATA.mkdir(parents=True, exist_ok=True)
    with CSV.open("w", encoding="utf-8", newline="") as f:
        f.write("n,lambda_n_trunc,tail_bound,lambda_over_half_n_log_n\n")
        for i in range(N_MAX):
            r = ratio[i] if n[i] >= 2 else ""
            f.write(f"{n[i]},{lam[i]:.10g},{tail[i]:.6g},{r}\n")
    meta = {
        "experiment": "E-008",
        "task": "3.2 part 1 (Li coefficients from zero table)",
        "description": "Li coefficients lambda_n (n<=2000) from E-002 zeros "
                       "with explicit tail bound; validation vs lambda_1 "
                       "literature value; leading-asymptotic ratio trend",
        "script": "research/numerics/scripts/E-008-li-coefficients.py",
        "parameters": {"n_max": N_MAX, "zeros_used": int(len(gam)),
                       "H": H, "precision": "float64"},
        "checks": checks,
        "key_values": {
            "lambda_1": float(lam[0]),
            "lambda_100": float(lam[99]),
            "lambda_1000": float(lam[999]),
            "lambda_2000": float(lam[1999]),
        },
        "usable_range_n": n_usable,
        "caveats": ["truncated at H; tail bound explicit in CSV; usable "
                    "range (tail<5%) is only n<=~116 with 10^4 zeros — "
                    "n~10^5 requires the prime-side Bombieri-Lagarias "
                    "exact formulas (noted as task 3.2 continuation)",
                    "positivity on our data is trivial (zeros on-line) — "
                    "not evidence about RH"],
        "elapsed_seconds": round(time.time() - t0, 1),
        "timestamp_utc": datetime.now(timezone.utc).isoformat(),
        "row_count": N_MAX,
    }
    JSON.write_text(json.dumps(meta, indent=2), encoding="utf-8")
    ok = all(checks.values())
    print(f"wrote {CSV.name}, {JSON.name} — "
          f"{'ALL CHECKS PASS' if ok else 'CHECK FAILURE'}")
    return 0 if ok else 1


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