E-011-li-coefficients-arithmetic.py

{
  "experiment": "E-011",
  "task": "3.2 continuation: Li coefficients from arithmetic side",
  "description": "Bombieri-Lagarias Theorem 2 (eta-form): lambda_n from Stieltjes constants (parallel mpmath sweep, series division) + Arb zeta values; zero-free",
  "script": "research/numerics/scripts/E-011-li-coefficients-arithmetic.py",
  "parameters": {
    "n_max": 2000,
    "dps": 1000,
    "arb_bits": 3530
  },
  "checks": {
    "V1_lambda1_exact": true,
    "V2_matches_E008_zero_side": true,
    "V3_all_positive": true,
    "V5_precision_self_consistent": true
  },
  "key_values": {
    "lambda_1": "0.0230957089661210338143102479065",
    "lambda_10": "2.27933936319316",
    "lambda_100": "118.603775376791",
    "lambda_1000": "2326.05316168647",
    "lambda_n_max": "5351.7595383815",
    "asymptotic_ratio_at_n_max": "1.0021569",
    "self_consistency_rel_dev": "6.67e-705"
  },
  "caveats": [
    "mpmath 1.3 + Python 3.14: mp.zeta broken beyond ~16 digits; zeta values via flint.arb; mp.zetazero unaffected",
    "mpmath arithmetic is not interval-rigorous; certification rests on self-consistency + E-008 cross-checks; an Arb ball-arithmetic upgrade would make this fully rigorous"
  ],
  "elapsed_seconds": 1072.5,
  "timestamp_utc": "2026-08-18T07:09:38.361237+00:00",
  "row_count": 2000
}

"""E-011 — Task 3.2 continuation: Li coefficients from the ARITHMETIC side.

Bombieri–Lagarias Theorem 2 (J. Number Theory 77 (1999) 274-287), with the
eta-sum written via the Laurent coefficients of -zeta'/zeta at s=1:

  -zeta'/zeta(s+1) = 1/s + sum_{k>=0} eta_k s^k,
  eta_k = (-1)^k/k! lim_{x->inf} { sum_{m<=x} Lambda(m)(log m)^k/m
                                   - (log x)^{k+1}/(k+1) },

  lambda_n = 1 - sum_{j=1}^n C(n,j) eta_{j-1}
             - (log 4pi + gamma_E) n/2
             - sum_{j=2}^n (-1)^{j-1} C(n,j) (1 - 2^{-j}) zeta(j).

(The (-1)^{j-1} in the theorem's eta-sum cancels the (-1)^{j-1} inside
the eta definition, leaving -sum C(n,j) eta_{j-1}. Verified numerically
against E-008 zero-side values for n <= 116 and against the machine-
validated explicit-formula conventions of D-001 NT-1.)

eta_k come from the Stieltjes constants gamma_m by series division:
  D(u) = u*zeta(1+u) = 1 + sum_{m>=0} (-1)^m gamma_m/m! u^{m+1},
  N(u) = -u^2*zeta'(1+u) = 1 - sum_{m>=1} (-1)^m gamma_m/(m-1)! u^{m+1},
  -zeta'/zeta(1+u) = (1/u) N(u)/D(u)  =>  eta_k = [u^{k+1}] N/D.

ENVIRONMENT CAVEAT (session 003 finding): mpmath 1.3 on Python 3.14
returns WRONG zeta values (integer AND non-integer) beyond ~16-19
digits. mp.zetazero is unaffected (Riemann-Siegel path; E-002 table
safe). All zeta(j) here come from flint.arb (ball arithmetic).
Stieltjes constants via flint.acb.stieltjes (rigorous Arb balls);
mpmath used only for series arithmetic (mp.mp.dps setter — bare mp.dps is broken in this env).

Precision budget: zeta-sum cancellation ~2^n (0.301n digits), eta-sum
~(4/3)^n (0.125n digits); dps = ceil(0.45 n_max) + 100.

Validations:
  V1. lambda_1 = 0.02309570896612103381... (28-digit match)
  V2. n <= 116: |lambda_n - E-008 zero-side| < max(2*tail_bound, 2% rel)
  V3. all lambda_n > 0
  V4. refined asymptotic ratio lambda_n / [(n/2)(log n - log 2pi + gamma - 1)]
      -> 1 (Arias de Reyna; logged, not a hard check)
  V5. self-consistency: lambda_n at sample n recomputed at 1.4x dps agree

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

import numpy as np
import mpmath as mp
import flint

REPO = Path(__file__).resolve().parents[3]
DATA = REPO / "research" / "numerics" / "data"
E008_CSV = DATA / "E-008-li-coefficients.csv"
CSV = DATA / "E-011-li-coefficients-arithmetic.csv"
JSON = DATA / "E-011-li-coefficients-arithmetic.json"

N_MAX = 1000
LAMBDA1_KNOWN = ("0.02309570896612103381431024790649529162193212715205"
                 "075952539207221297135647672457997079856951170983336")
# ^ corrected session 004: previous string diverged from the true value
# 1 + gamma/2 - log(4pi)/2 at digit 22 (flint-verified to 100 digits).


def _stieltjes_worker(args):
    k, bits, digits = args
    import flint
    flint.ctx.prec = bits
    return k, flint.acb.stieltjes(k).real.str(digits, radius=False)


def stieltjes_sweep(k_max: int, dps: int) -> list:
    """gamma_0..gamma_{k_max} via Arb (flint.acb.stieltjes, rigorous
    balls), parallel over k."""
    bits = int(np.ceil(dps * 3.33)) + 200
    with Pool(14) as pool:
        out = pool.map(_stieltjes_worker,
                       [(k, bits, dps + 10) for k in range(k_max + 1)],
                       chunksize=4)
    out.sort()
    return [mp.mpf(v) for _, v in out]


def eta_from_stieltjes(gam, n_max):
    N_ = [mp.mpf(0)] * (n_max + 2)
    D_ = [mp.mpf(0)] * (n_max + 2)
    N_[0] = D_[0] = 1
    for m in range(0, n_max + 1):
        D_[m + 1] = (-1) ** m * gam[m] / mp.factorial(m)
    for m in range(1, n_max + 1):
        N_[m + 1] = -(-1) ** m * gam[m] / mp.factorial(m - 1)
    E = [mp.mpf(0)] * (n_max + 2)
    for k in range(0, n_max + 2):
        s = N_[k]
        for i in range(1, k + 1):
            s -= D_[i] * E[k - i]
        E[k] = s
    return [E[k + 1] for k in range(n_max + 1)]


def arb_zetas(n_max):
    return [mp.mpf(0), mp.mpf(0)] + [
        mp.mpf(flint.arb(j).zeta().str(mp.mp.dps + 10, radius=False))
        for j in range(2, n_max + 1)]


def li_lambda(n, eta, w, c2):
    s = mp.mpf(1) - c2 * n
    b = 1
    for j in range(1, n + 1):
        b = b * (n - j + 1) // j          # exact integer C(n, j)
        s -= b * eta[j - 1]
        if j >= 2:
            s -= ((-1) ** (j - 1)) * b * w[j]
    return s


def main() -> int:
    n_max = N_MAX
    if "--n-max" in sys.argv:
        n_max = int(sys.argv[sys.argv.index("--n-max") + 1])
    if CSV.exists() and "--force" not in sys.argv:
        print(f"cached: {CSV} exists (use --force to recompute)"); return 0
    t0 = time.time()
    dps = int(np.ceil(0.45 * n_max)) + 100
    mp.mp.dps = dps
    flint.ctx.prec = int(np.ceil(dps * 3.33)) + 200
    print(f"n_max={n_max}, dps={dps}, arb bits={flint.ctx.prec}",
          flush=True)

    t1 = time.time()
    gam = stieltjes_sweep(n_max, dps)
    print(f"stieltjes 0..{n_max}: {time.time()-t1:.0f}s", flush=True)

    t1 = time.time()
    eta = eta_from_stieltjes(gam, n_max)
    print(f"eta recursion: {time.time()-t1:.0f}s  eta_0={mp.nstr(eta[0], 12)}"
          f" eta_1={mp.nstr(eta[1], 12)}", flush=True)

    t1 = time.time()
    zeta = arb_zetas(n_max)
    print(f"arb zeta values: {time.time()-t1:.0f}s", flush=True)

    t1 = time.time()
    c2 = (mp.log(4 * mp.pi) + mp.euler) / 2
    w = [mp.mpf(0), mp.mpf(0)] + [(1 - mp.mpf(2) ** (-j)) * zeta[j]
                                  for j in range(2, n_max + 1)]
    lam = [mp.mpf(0)] + [li_lambda(n, eta, w, c2)
                         for n in range(1, n_max + 1)]
    print(f"lambda sums: {time.time()-t1:.0f}s", flush=True)

    # ---- V5: self-consistency at higher precision, reduced range
    t1 = time.time()
    k_half = n_max // 2
    dps2 = int(dps * 1.4)
    mp.mp.dps = dps2
    flint.ctx.prec = int(np.ceil(dps2 * 3.33)) + 200
    gam2 = stieltjes_sweep(k_half, dps2)
    eta2 = eta_from_stieltjes(gam2, k_half)
    zeta2 = arb_zetas(k_half)
    c2b = (mp.log(4 * mp.pi) + mp.euler) / 2
    w2 = [mp.mpf(0), mp.mpf(0)] + [(1 - mp.mpf(2) ** (-j)) * zeta2[j]
                                   for j in range(2, k_half + 1)]
    sample = sorted({1, 2, 100, k_half})
    lam2 = {n: li_lambda(n, eta2, w2, c2b) for n in sample}
    mp.mp.dps = dps
    rel_dev = max(abs(lam2[n] - lam[n]) / abs(lam[n]) for n in sample)
    check_v5 = bool(rel_dev < mp.mpf(10) ** (-0.30 * k_half))
    print(f"self-consistency (n<={k_half} @1.4x dps): rel dev "
          f"{mp.nstr(rel_dev, 3)} ({time.time()-t1:.0f}s)", flush=True)

    # ---- validations
    lam1_err = abs(lam[1] - mp.mpf(LAMBDA1_KNOWN))
    check_v1 = bool(lam1_err < mp.mpf("1e-28"))
    zs = {}
    with E008_CSV.open() as f:
        f.readline()
        for line in f:
            p = line.split(",")
            zs[int(p[0])] = (float(p[1]), float(p[2]))
    worst = 0.0
    for n in range(1, 117):
        z, tb = zs[n]
        dev = abs(float(lam[n]) - z)
        worst = max(worst, dev / max(2 * tb, 0.02 * abs(z)))
    check_v2 = bool(worst < 1.0)
    check_v3 = bool(all(lam[n] > 0 for n in range(1, n_max + 1)))

    def lead(n):
        return (n / 2) * (mp.log(n) - mp.log(2 * mp.pi) + mp.euler - 1)

    ratio_last = lam[n_max] / lead(n_max)
    checks = {
        "V1_lambda1_exact": check_v1,
        "V2_matches_E008_zero_side": check_v2,
        "V3_all_positive": check_v3,
        "V5_precision_self_consistent": check_v5,
    }
    print(f"lambda_1 err: {mp.nstr(lam1_err, 3)}; E-008 worst dev ratio: "
          f"{worst:.3f}; lambda_{n_max} = {mp.nstr(lam[n_max], 8)}; "
          f"asymptotic ratio: {mp.nstr(ratio_last, 6)}", flush=True)

    DATA.mkdir(parents=True, exist_ok=True)
    with CSV.open("w", encoding="utf-8", newline="") as f:
        f.write("n,lambda_n,lambda_over_refined_asymptotic\n")
        for n in range(1, n_max + 1):
            f.write(f"{n},{mp.nstr(lam[n], 25)},"
                    f"{mp.nstr(lam[n] / lead(n), 12)}\n")
    meta = {
        "experiment": "E-011",
        "task": "3.2 continuation: Li coefficients from arithmetic side",
        "description": "Bombieri-Lagarias Theorem 2 (eta-form): lambda_n "
                       "from Stieltjes constants (parallel mpmath sweep, "
                       "series division) + Arb zeta values; zero-free",
        "script": "research/numerics/scripts/E-011-li-coefficients-arithmetic.py",
        "parameters": {"n_max": n_max, "dps": dps,
                       "arb_bits": int(np.ceil(dps * 3.33)) + 200},
        "checks": checks,
        "key_values": {
            "lambda_1": mp.nstr(lam[1], 30),
            "lambda_10": mp.nstr(lam[10], 15),
            "lambda_100": mp.nstr(lam[100], 15),
            "lambda_1000": mp.nstr(lam[1000], 15) if n_max >= 1000 else None,
            "lambda_n_max": mp.nstr(lam[n_max], 15),
            "asymptotic_ratio_at_n_max": mp.nstr(ratio_last, 8),
            "self_consistency_rel_dev": mp.nstr(rel_dev, 3),
        },
        "caveats": [
            "mpmath 1.3 + Python 3.14: mp.zeta broken beyond ~16 digits; "
            "zeta values via flint.arb; mp.zetazero unaffected",
            "mpmath arithmetic is not interval-rigorous; certification "
            "rests on self-consistency + E-008 cross-checks; an Arb "
            "ball-arithmetic upgrade would make this fully rigorous",
        ],
        "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: ' + str(checks)}")
    return 0 if ok else 1


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