E-028-traced-margin.py

{
  "experiment": "E-028",
  "task": "4.5a support: D-056 traced R-bound verification + margin rerun with A=10.5, B0=48 (honest C0 incl. log log term)",
  "script": "research/numerics/scripts/E-028-traced-margin.py",
  "checks": {
    "V_R_traced_bound_holds": {
      "pass": true,
      "points": [
        {
          "sigma": -2.5,
          "t": 10.0,
          "abs_R": 0.5364627998701963,
          "bound": 74.091519822774,
          "holds": true
        },
        {
          "sigma": -2.5,
          "t": 50.0,
          "abs_R": 1.8035543334872959,
          "bound": 89.48805904510499,
          "holds": true
        },
        {
          "sigma": -2.5,
          "t": 100.0,
          "abs_R": 2.726827160278793,
          "bound": 96.56221453948484,
          "holds": true
        },
        {
          "sigma": -2.5,
          "t": 500.0,
          "abs_R": 4.0305314954089475,
          "bound": 113.29530125676315,
          "holds": true
        },
        {
          "sigma": -2.5,
          "t": 1000.0,
          "abs_R": 4.7507714994958565,
          "bound": 120.5524094572705,
          "holds": true
        },
        {
          "sigma": 0.5,
          "t": 10.0,
          "abs_R": 0.232976073221556,
          "bound": 74.091519822774,
          "holds": true
        },
        {
          "sigma": 0.5,
          "t": 50.0,
          "abs_R": 1.0575019981720868,
          "bound": 89.48805904510499,
          "holds": true
        },
        {
          "sigma": 0.5,
          "t": 100.0,
          "abs_R": 1.3861083473172386,
          "bound": 96.56221453948484,
          "holds": true
        },
        {
          "sigma": 0.5,
          "t": 500.0,
          "abs_R": 2.191750394872167,
          "bound": 113.29530125676315,
          "holds": true
        },
        {
          "sigma": 0.5,
          "t": 1000.0,
          "abs_R": 2.535043296015203,
          "bound": 120.5524094572705,
          "holds": true
        },
        {
          "sigma": 2.0,
          "t": 10.0,
          "abs_R": 0.11614405022018243,
          "bound": 74.091519822774,
          "holds": true
        },
        {
          "sigma": 2.0,
          "t": 50.0,
          "abs_R": 0.47282316679931,
          "bound": 89.48805904510499,
          "holds": true
        },
        {
          "sigma": 2.0,
          "t": 100.0,
          "abs_R": 0.13345107201883163,
          "bound": 96.56221453948484,
          "holds": true
        },
        {
          "sigma": 2.0,
          "t": 500.0,
          "abs_R": 0.6914268199229783,
          "bound": 113.29530125676315,
          "holds": true
        },
        {
          "sigma": 2.0,
          "t": 1000.0,
          "abs_R": 0.6591087549534896,
          "bound": 120.5524094572705,
          "holds": true
        }
      ]
    },
    "V5_Z2_certified": {
      "pass": true,
      "partial": 0.5699599927023044,
      "tail": 1.4815510557964273e-05,
      "bound": 0.5699748082128623
    },
    "V2_margin_positive_convex": true
  },
  "key_values": {
    "M0_from_E023": 1.335180201468999,
    "N0_traced_10p5_48": 75451,
    "margin_N0": 2.396310426061973,
    "overlap_factor_vs_MV_5e5": 6.626817404673232,
    "margin_at": {
      "75451": 2.396310426061973,
      "100000": 92785.67601524561,
      "500000": 2006008.571308432,
      "1000000": 4755774.032035855,
      "1000000000": 8921657382.546026,
      "1000000000000": 12395765983564.96,
      "1000000000000000000": 1.9304326196882584e+19,
      "8999999999994000000000001": 2.4579625383799396e+26
    },
    "N1": 8999999999994000000000001
  },
  "elapsed_seconds": 3.7
}

"""E-028 — task 4.5a support: D-056 traced R-bound verification + margin rerun.

Components:
  V-R. Verify the traced reference-point R-bound (D-056 Theorem R)
       |R(s)| <= 6.0 log(|t|+2) + 40,  sigma in [-5/2, 3], |t| >= 10
       at sample points sigma in {-2.5, 0.5, 2}, t in {10,50,100,500,1000}
       (R(s) = zeta'/zeta(s) - sum_{|gamma-t|<1} 1/(s-rho), zeta'/zeta
       via the E-021 xi-route at dps 30; zeros from the E-027 cache).
  V5.  Certified Z2 = -zeta'/zeta(2) <= 0.5699: partial sum to 1e4 plus
       tail bound (log N + 1)/N.
  S3'. Margin rerun with the TRACED constants A = 6.0, B0 = 40 (was:
       measured 0.9/1.0, first-estimate 2.0/2.0) and the I_I length fix
       2(sqrt(n)+1). Recompute N0, contiguity to 4e5, convexity grid to
       N1, margin references. M0 = 1.3352 imported from E-023 (S1).

Outputs: research/numerics/data/E-028-traced-margin.csv (+ .json)
"""
import csv as csvmod
import json
import sys
import time
from pathlib import Path

import mpmath as mp
import numpy as np

REPO = Path(__file__).resolve().parents[3]
DATA = REPO / "research" / "numerics" / "data"
E021 = DATA / "E-021-sinf-sf-curve.csv"
E023J = DATA / "E-023-assembly.json"
ZCACHE = DATA / "E-027-zeros-cache.npy"
CSV = DATA / "E-028-traced-margin.csv"
JSON = DATA / "E-028-traced-margin.json"

GAMMA = 0.577215664901532860606512090082402431
C1L = (GAMMA - 1 - np.log(2 * np.pi)) / 2 - GAMMA / 2
A_TR, B0_TR = 10.5, 48.0
T0_PLATT = 3.0e12


def zeta_logder(s):
    sm = mp.mpc(s.real, s.imag)

    def xi_mp(u):
        return 0.5 * u * (u - 1) * mp.pi ** (-u / 2) * mp.gamma(u / 2) \
            * mp.zeta(u)

    xi_log = mp.diff(xi_mp, sm) / xi_mp(sm)
    return complex(xi_log - 1 / sm - 1 / (sm - 1)
                   + mp.mpf(0.5) * mp.log(mp.pi)
                   - mp.mpf(0.5) * mp.digamma(sm / 2))


def check_R_bound(zeros):
    mp.mp.dps = 30
    out = []
    for sigma in (-2.5, 0.5, 2.0):
        for t in (10.0, 50.0, 100.0, 500.0, 1000.0):
            s = complex(sigma, t)
            val = zeta_logder(s)
            near = [g for g in zeros if abs(g - t) < 1.0]
            r = val - sum(1.0 / (s - complex(0.5, g)) for g in near)
            bound = 10.5 * np.log(t + 2) + 48.0
            out.append({"sigma": sigma, "t": t, "abs_R": abs(r),
                        "bound": bound, "holds": bool(abs(r) <= bound)})
    return out


def certified_Z2(N=10 ** 6):
    """-zeta'/zeta(2) = sum Lambda(n)/n^2; partial sum + tail bound
    (log N + 1)/N. True value 0.569955...; certify <= 0.5700."""
    lam = np.zeros(N + 1)
    sieve = np.zeros(N + 1, dtype=np.int64)
    for p in range(2, N + 1):
        if sieve[p] == 0:
            sieve[p::p] = p
    for n in range(2, N + 1):
        p = sieve[n]
        if p > 0 and abs(p ** round(np.log(n) / np.log(p)) - n) < 0.5:
            lam[n] = np.log(p)
    part = float(np.sum(lam[2:] / np.arange(2, N + 1) ** 2))
    tail = (np.log(N) + 1) / N
    return part, tail, part + tail


def F_count(T):
    return (T / (2 * np.pi)) * np.log(T / (2 * np.pi)) - T / (2 * np.pi) \
        + 7 / 8


def S_bound(T):
    return 0.137 * np.log(T) + 0.443 * np.log(np.log(T)) + 1.588


def C0(t):
    return (F_count(t + 1) - F_count(t)) + S_bound(t + 1) + S_bound(t)


def H_bound(T):
    return 3 * C0(T) ** 2 + A_TR * np.log(T + 2) + B0_TR


def B_bound(n, m0):
    n = np.asarray(n, dtype=np.float64)
    T = np.sqrt(n) + 1.0
    twoN = 2 * (F_count(T) + S_bound(T))
    I3_long = 4 * T * ((2 / 3) * C0(T) + A_TR * np.log(T + 2) + B0_TR)
    I3_short = 40.0 * m0
    I24 = 2 * (np.e ** 4 + 1) * (7 * H_bound(T) + 3.3)
    r = np.sqrt(n)
    # I_I length fix: 2T = 2(sqrt(n)+1)
    logI1 = np.log(2 * (r + 1) * 0.5700) + r / 2 - (2 * r - 2) * np.log(2)
    I1 = np.where(logI1 > -700, np.exp(np.maximum(logI1, -745)), 0.0)
    return twoN + I3_long + I3_short + I24 + I1 + 1.0


def margin(n, m0):
    n = np.asarray(n, dtype=np.float64)
    return 0.5 * n * np.log(n) + C1L * n - 4.0 - B_bound(n, m0) + 1.0


def main() -> int:
    if JSON.exists() and "--force" not in sys.argv:
        print(f"cached: {JSON} (use --force to recompute)")
        return 0
    t0 = time.time()

    m0 = json.loads(E023J.read_text())["key_values"]["M0_use"]
    print(f"M0 imported from E-023: {m0}")

    print("V-R: traced R-bound verification ...", flush=True)
    zeros = np.load(ZCACHE)
    rpts = check_R_bound(zeros)
    vr = all(p["holds"] for p in rpts)
    worst = max(rpts, key=lambda p: p["abs_R"] / p["bound"])
    print(f"  all {len(rpts)} points hold: {vr}; worst ratio "
          f"{worst['abs_R']/worst['bound']:.3f} at "
          f"sigma={worst['sigma']}, t={worst['t']}", flush=True)

    print("V5: certified Z2 ...", flush=True)
    part, tail, z2 = certified_Z2()
    v5 = bool(z2 <= 0.5700)
    print(f"  partial {part:.6f} + tail {tail:.6f} = {z2:.6f} "
          f"<= 0.5700: {v5}", flush=True)

    print("S3': margin rerun with traced constants ...", flush=True)
    ngrid = np.arange(20, 400001)
    marg = margin(ngrid, m0)
    pos = ngrid[(marg > 0) & (np.gradient(marg, ngrid) > 0)]
    N0 = int(pos[0]) if len(pos) else -1
    cont = bool(np.all(marg[ngrid >= N0] > 0)) if N0 > 0 else False
    margN0 = float(marg[ngrid == N0][0]) if N0 > 0 else float("nan")
    print(f"  N0 (traced 6.0/40) = {N0} (contiguous to 4e5: {cont}, "
          f"margin(N0) = {margN0:.1f})", flush=True)

    logn = np.unique(np.logspace(np.log10(max(N0, 20)), np.log10(9e24),
                                 4000).astype(np.float64))
    bmarg = margin(logn, m0)
    conv = bool(np.all(np.diff(bmarg, 2) > -1e-6 * bmarg[:-2].max()))

    N1 = int(3_000_000_000_000 - 1) ** 2
    refs = {int(n): float(margin(n, m0)) for n in
            [N0, 100000, 500000, 10 ** 6, 10 ** 9, 10 ** 12, 10 ** 18,
             N1]}
    for k, v in refs.items():
        print(f"  margin({k:.3e}) = {v:.4e}", flush=True)
    v2 = bool(N0 > 0 and cont and conv and N0 < 500000)
    print(f"  convexity grid: {conv}; N0 < 5e5 (MV coverage): "
          f"{N0 < 500000}", flush=True)

    with open(CSV, "w", newline="") as f:
        w = csvmod.writer(f)
        w.writerow(["n", "margin_traced_6_40"])
        for n in range(20, 400001, 20):
            w.writerow([n, f"{float(margin(n, m0)):.6f}"])

    meta = {
        "experiment": "E-028",
        "task": "4.5a support: D-056 traced R-bound verification + "
                "margin rerun with A=10.5, B0=48 (honest C0 incl. log log term)",
        "script": "research/numerics/scripts/E-028-traced-margin.py",
        "checks": {
            "V_R_traced_bound_holds": {"pass": vr, "points": rpts},
            "V5_Z2_certified": {"pass": v5, "partial": part,
                                "tail": tail, "bound": z2},
            "V2_margin_positive_convex": v2,
        },
        "key_values": {
            "M0_from_E023": m0,
            "N0_traced_10p5_48": N0,
            "margin_N0": margN0,
            "overlap_factor_vs_MV_5e5": 500000 / N0 if N0 > 0 else None,
            "margin_at": {str(k): v for k, v in refs.items()},
            "N1": N1,
        },
        "elapsed_seconds": round(time.time() - t0, 1),
    }
    JSON.write_text(json.dumps(meta, indent=2))
    print(f"  wrote {CSV.name}, {JSON.name}; elapsed "
          f"{meta['elapsed_seconds']}s")
    return 0 if (vr and v5 and v2) else 1


if __name__ == "__main__":
    sys.exit(main())