E-033-repaired-stepping-validation.py

{
  "experiment": "E-033",
  "title": "F-D side thread 6-R-1: cross-validation of the repaired-stepping constants (D-073 secs.4-6): Z0 certified, q-table, R(c) vs E-032 measured, kernel charge, fold fit, margins at c = 1.5 and 2",
  "session": "039",
  "thread": "ii-prime-6-R",
  "design_doc": "D-073",
  "checks": {
    "V1_Z0_certified": true,
    "V2_q_table_and_boundary": true,
    "V3_R_covers_E032_measured": true,
    "V4_kernel_charge": true,
    "V5_fold_fit": true,
    "V6_margin_positive": true,
    "V7_full_fold": true
  },
  "all_checks_pass": true,
  "dps": 50,
  "constants": {
    "Z0_lo": "0.569921006573",
    "Z0_hi": "0.570366088021",
    "q_1.5": "0.77981175",
    "q_2": "0.64201271",
    "R_1.5": "5.9392431",
    "R_2": "1.9734828",
    "kernel_1.5": "6.9166936",
    "kernel_2": "3.7182818",
    "N1_1.5": "4.0e+24",
    "N1_2": "2.25e+24"
  },
  "notes": "Recommendation (D-073 sec.5): c = 2 for the main line (R = 1.97 < 3.3 broken charge, kernel 3.72 < 55.6, N1 = 2.25e24); c = 1.5 the range option (N1 = 4e24, R = 5.94, kernel 6.92)."
}

#!/usr/bin/env python3
"""E-033 — F-D side thread ii′-6-R, sub-stage 6-R-1 (deliverable D2-v):
cross-validation of ALL constants in the repaired horizontal proof
(D-073 §§4–6): the corrected stepping lemma at T = c.sqrt(n),
c > 1/sqrt(log 2), with the explicit right-tail bound

    R(c) := Z0 . [ e^{4/c^2} / (4(1-q) log 2) + 1/(2 log 2) ],
    q := e^{1/c^2}/2,  Z0 := |zeta'/zeta(2)|.

Checks:
  V1: Z0 certified: partial Dirichlet sum + certified tail bound
      bracketing |zeta'/zeta(2)| (mpmath cross-check).
  V2: the q-table: q(1) > 1 (the F-D bug), q(c) < 1 for c in
      {1.25, 1.3, 1.5, 1.75, 2}; boundary c* = 1/sqrt(log 2) in
      [1.2011, 1.2013].
  V3: R(c) covers the E-032 MEASURED |I_right| at c = 1.5, 2
      (reads the committed E-032 CSV; n = 1e3, 1e4, 1e5).
  V4: kernel charge: max_{sigma in [-3,3]} |k_n(sigma+iT)| <= e^{4/c^2}+1
      at n = 1e5, T = c.sqrt(n), c in {1.5, 2} (direct evaluation,
      0.02 grid).
  V5: horizontal fold fit at T = c.sqrt(n)+1:
      NI(T) + 5.5(10.5 log(T+2)+48) <= 7.Hhat(T) on the grid,
      c in {1.5, 2} (the E-030 V3 analogue at the repaired height).
  V6: margin(n; c) = sInftyLower(n) - Bn_c(n) + 1 > 0 on a dense log
      grid over [8e4, N1(c)], N1(c) = ((3e12-1)/c)^2, c in {1.5, 2}.
  V7: full fold honest_c(n) <= Bn_c(n) on the grid (E-030 V4 analogue).

All numerics dps 50, deterministic; reads E-032 CSV (committed data).
"""
import csv as csvmod
import json
import sys
from pathlib import Path

import mpmath as mp

mp.mp.dps = 50

DATA = Path(__file__).resolve().parent.parent / "data"
CSV = DATA / "E-033-repaired-stepping-validation.csv"
JSON = DATA / "E-033-repaired-stepping-validation.json"
E032CSV = DATA / "E-032-horizontal-true-integral.csv"

pi, e = mp.pi, mp.e
LOG2 = mp.log(2)

SIEVE_MAX = 25000
_spf = list(range(SIEVE_MAX + 1))
for _i in range(2, int(SIEVE_MAX ** 0.5) + 1):
    if _spf[_i] == _i:
        for _j in range(_i * _i, SIEVE_MAX + 1, _i):
            if _spf[_j] == _j:
                _spf[_j] = _i
LAM = [mp.mpf(0)] * (SIEVE_MAX + 1)
LOGK = [mp.mpf(0)] * (SIEVE_MAX + 1)
for _k in range(2, SIEVE_MAX + 1):
    LOGK[_k] = mp.log(_k)
    _p = _spf[_k]
    _m = _k
    while _m % _p == 0:
        _m //= _p
    if _m == 1:
        LAM[_k] = LOGK[_p]


def Z0_bracket(K=25000):
    """|zeta'/zeta(2)| = sum Lambda(k)/k^2, certified tail
    <= integral_{K-1}^inf log x / x^2 dx = (log(K-1)+1)/(K-1)."""
    s = mp.mpf(0)
    for k in range(2, K + 1):
        if LAM[k]:
            s += LAM[k] / (mp.mpf(k) ** 2)
    tail = (mp.log(K - 1) + 1) / (K - 1)
    return s, s + tail


def qc(c):
    return e ** (1 / mp.mpf(c) ** 2) / 2


def Rc(c, Z0hi):
    c = mp.mpf(c)
    q = qc(c)
    return Z0hi * (e ** (4 / c ** 2) / (4 * (1 - q) * LOG2) + 1 / (2 * LOG2))


def kernel_charge(c):
    return e ** (4 / mp.mpf(c) ** 2) + 1


def C0hat(t):
    return mp.mpf("0.433") * mp.log(t + 2) + mp.mpf("0.886") * mp.log(mp.log(t + 2)) + mp.mpf("2.9")


def Z2(t):
    return mp.mpf("0.6") * mp.log(t + 2) + mp.mpf("0.886") * mp.log(mp.log(t + 2)) + mp.mpf("2.7")


def FAlt(T):
    return ((1 / (2 * pi)) * (T * mp.log(T)) + mp.mpf("0.137") * mp.log(T)
            + mp.mpf("0.443") * mp.log(mp.log(T)) + (mp.mpf(7) / 8 + mp.mpf("1.588")))


def Hhat(T):
    return 3 * C0hat(T) ** 2 + mp.mpf("10.5") * mp.log(T + 2) + 48


def NI(T, a):
    """Near-sum integral on the horizontal at height T, pigeonhole
    window [a, a+1], clearance 1/(9.C0hat(a+2))."""
    c_clear = 1 / (9 * C0hat(a + 2))
    return Z2(T) * 2 * mp.log(mp.mpf("16.5") / c_clear)


def IIMaj_c(n, c):
    return 2 * (c * mp.sqrt(n) + 1) * (mp.e ** (mp.sqrt(n) / 2) + 1) * mp.mpf("0.57") * 2 ** (-(2 * mp.sqrt(n) - 2))


def Bn_c(n, c, Z0hi):
    T = c * mp.sqrt(n) + 1
    return (4 * FAlt(T) + 4 * T * ((mp.mpf(2) / 3) * Z2(T) + mp.mpf("10.5") * mp.log(T + 2) + 48)
            + 40 * mp.mpf("1.3352")
            + 2 * kernel_charge(c) * (7 * Hhat(T) + Rc(c, Z0hi))
            + IIMaj_c(n, c) + 1)


def honest_c(n, c, Z0hi):
    T = c * mp.sqrt(n) + 1
    a = c * mp.sqrt(n)
    horiz = 2 * kernel_charge(c) * (NI(T, a) + mp.mpf("5.5") * (mp.mpf("10.5") * mp.log(T + 2) + 48)
                                    + Rc(c, Z0hi))
    return (4 * FAlt(T) + 4 * T * ((mp.mpf(2) / 3) * Z2(T) + mp.mpf("10.5") * mp.log(T + 2) + 48)
            + 40 * mp.mpf("1.3352") + horiz + IIMaj_c(n, c) + 1)


def sInftyLower(n):
    return mp.mpf("0.5") * n * mp.log(n) + (mp.mpf("-1.41894")) * n - 4


def N1(c):
    return ((mp.mpf("3e12") - 1) / c) ** 2


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

    rows = []

    # V1: Z0 certified
    zlo, zhi = Z0_bracket()
    zmp = -mp.zeta(2, derivative=1) / mp.zeta(2)
    v1 = zlo <= zmp <= zhi and mp.mpf("0.5695") < zlo and zhi < mp.mpf("0.5705") \
        and zhi - zlo < mp.mpf("1e-3")
    rows.append({"check": "V1_Z0_certified", "c": "", "n": "",
                 "value_lo": mp.nstr(zlo, 12), "value_hi": mp.nstr(zhi, 12),
                 "mpmath": mp.nstr(zmp, 12), "pass": str(bool(v1))})

    # V2: q-table + boundary
    cs = ["1", "1.25", "1.3", "1.5", "1.75", "2"]
    qtab = {c: qc(mp.mpf(c)) for c in cs}
    cstar = 1 / mp.sqrt(LOG2)
    v2 = qtab["1"] > 1 and all(qtab[c] < 1 for c in cs[1:]) \
        and mp.mpf("1.2011") < cstar < mp.mpf("1.2013")
    for c in cs:
        rows.append({"check": "V2_q_table", "c": c, "n": "",
                     "value_lo": mp.nstr(qtab[c], 6),
                     "value_hi": mp.nstr(Rc(mp.mpf(c), zhi), 6),
                     "mpmath": mp.nstr(kernel_charge(mp.mpf(c)), 6) if c != "1" else "",
                     "pass": str(bool(qtab[c] < 1))})
    rows.append({"check": "V2_c_boundary", "c": "c*", "n": "",
                 "value_lo": mp.nstr(cstar, 8), "value_hi": "1/sqrt(log2)",
                 "mpmath": "", "pass": str(bool(mp.mpf("1.2011") < cstar < mp.mpf("1.2013")))})

    # V3: R(c) covers the E-032 measured |I_right|
    measured = {}
    with E032CSV.open() as f:
        for r in csvmod.DictReader(f):
            if r["check"] == "right_part" and r["c"] in ("1.5", "2.0"):
                measured[(int(r["n"]), mp.mpf(r["c"]))] = mp.mpf(r["abs_I_right"])
    v3 = True
    for (n, c), val in sorted(measured.items()):
        ok = Rc(c, zhi) >= val
        v3 = v3 and ok
        rows.append({"check": "V3_R_covers_measured", "c": mp.nstr(c, 3), "n": str(n),
                     "value_lo": mp.nstr(val, 6), "value_hi": mp.nstr(Rc(c, zhi), 6),
                     "mpmath": "", "pass": str(bool(ok))})

    # V4: kernel charge on [-3,3]
    n = 10 ** 5
    v4 = True
    for c in [mp.mpf("1.5"), mp.mpf("2")]:
        T = c * mp.sqrt(n)
        mx = mp.mpf(0)
        sig = mp.mpf(-3)
        while sig <= 3:
            v = abs((1 + 1 / mp.mpc(sig, T)) ** n - 1)
            if v > mx:
                mx = v
            sig += mp.mpf("0.02")
        ok = mx <= kernel_charge(c)
        v4 = v4 and ok
        rows.append({"check": "V4_kernel_charge", "c": mp.nstr(c, 3), "n": str(n),
                     "value_lo": mp.nstr(mx, 6), "value_hi": mp.nstr(kernel_charge(c), 6),
                     "mpmath": "", "pass": str(bool(ok))})

    # grids
    def loggrid(lo, hi, pts):
        llo, lhi = mp.log(lo), mp.log(hi)
        return [mp.e ** (llo + (lhi - llo) * i / (pts - 1)) for i in range(pts)]

    # V5/V6/V7 per c
    v5, v6, v7 = True, True, True
    for c in [mp.mpf("1.5"), mp.mpf("2")]:
        grid = [mp.mpf("8e4"), mp.mpf("1e5")] + loggrid(mp.mpf("1e5"), N1(c), 30)
        maxfit = mp.mpf(0)
        minmargin = None
        for nn in grid:
            T = c * mp.sqrt(nn) + 1
            a = c * mp.sqrt(nn)
            fit = (NI(T, a) + mp.mpf("5.5") * (mp.mpf("10.5") * mp.log(T + 2) + 48)) / (7 * Hhat(T))
            if fit > maxfit:
                maxfit = fit
            v5 = v5 and (fit < 1)
            marg = sInftyLower(nn) - Bn_c(nn, c, zhi) + 1
            if minmargin is None or marg < minmargin:
                minmargin = marg
            v6 = v6 and (marg > 0)
            v7 = v7 and (honest_c(nn, c, zhi) <= Bn_c(nn, c, zhi))
        rows.append({"check": "V5_fold_fit_max_ratio", "c": mp.nstr(c, 3), "n": "grid",
                     "value_lo": mp.nstr(maxfit, 4), "value_hi": "< 1",
                     "mpmath": "", "pass": str(bool(maxfit < 1))})
        rows.append({"check": "V6_margin_min_on_grid", "c": mp.nstr(c, 3), "n": "grid",
                     "value_lo": mp.nstr(minmargin, 6), "value_hi": "> 0",
                     "mpmath": "", "pass": str(bool(minmargin > 0))})
        for nn in [mp.mpf("8e4"), mp.mpf("1e5"), N1(c)]:
            marg = sInftyLower(nn) - Bn_c(nn, c, zhi) + 1
            rows.append({"check": "V6_margin_value", "c": mp.nstr(c, 3), "n": mp.nstr(nn, 8),
                         "value_lo": mp.nstr(marg, 8), "value_hi": "> 0",
                         "mpmath": "", "pass": str(bool(marg > 0))})
        rows.append({"check": "V7_full_fold", "c": mp.nstr(c, 3), "n": "grid",
                     "value_lo": "honest<=Bn", "value_hi": "",
                     "mpmath": "", "pass": "see above"})

    # recompute v7 cleanly (the row above is a placeholder)
    v7 = all(honest_c(nn, c, zhi) <= Bn_c(nn, c, zhi)
             for c in [mp.mpf("1.5"), mp.mpf("2")]
             for nn in [mp.mpf("8e4"), mp.mpf("1e5")] + loggrid(mp.mpf("1e5"), N1(c), 30))

    allpass = v1 and v2 and v3 and v4 and v5 and v6 and v7

    DATA.mkdir(exist_ok=True)
    with CSV.open("w") as f:
        f.write(",".join(rows[0].keys()) + "\n")
        for r in rows:
            f.write(",".join(r.values()) + "\n")
    sidecar = {
        "experiment": "E-033",
        "title": "F-D side thread 6-R-1: cross-validation of the repaired-stepping "
                 "constants (D-073 secs.4-6): Z0 certified, q-table, R(c) vs E-032 "
                 "measured, kernel charge, fold fit, margins at c = 1.5 and 2",
        "session": "039",
        "thread": "ii-prime-6-R",
        "design_doc": "D-073",
        "checks": {"V1_Z0_certified": bool(v1), "V2_q_table_and_boundary": bool(v2),
                   "V3_R_covers_E032_measured": bool(v3), "V4_kernel_charge": bool(v4),
                   "V5_fold_fit": bool(v5), "V6_margin_positive": bool(v6),
                   "V7_full_fold": bool(v7)},
        "all_checks_pass": bool(allpass),
        "dps": 50,
        "constants": {"Z0_lo": mp.nstr(zlo, 12), "Z0_hi": mp.nstr(zhi, 12),
                      "q_1.5": mp.nstr(qc(mp.mpf("1.5")), 8),
                      "q_2": mp.nstr(qc(mp.mpf("2")), 8),
                      "R_1.5": mp.nstr(Rc(mp.mpf("1.5"), zhi), 8),
                      "R_2": mp.nstr(Rc(mp.mpf("2"), zhi), 8),
                      "kernel_1.5": mp.nstr(kernel_charge(mp.mpf("1.5")), 8),
                      "kernel_2": mp.nstr(kernel_charge(mp.mpf("2")), 8),
                      "N1_1.5": mp.nstr(N1(mp.mpf("1.5")), 6),
                      "N1_2": mp.nstr(N1(mp.mpf("2")), 6)},
        "notes": "Recommendation (D-073 sec.5): c = 2 for the main line "
                 "(R = 1.97 < 3.3 broken charge, kernel 3.72 < 55.6, N1 = 2.25e24); "
                 "c = 1.5 the range option (N1 = 4e24, R = 5.94, kernel 6.92).",
    }
    with JSON.open("w") as f:
        json.dump(sidecar, f, indent=2)
    print(f"V1 {bool(v1)} V2 {bool(v2)} V3 {bool(v3)} V4 {bool(v4)} "
          f"V5 {bool(v5)} V6 {bool(v6)} V7 {bool(v7)}")
    print(f"wrote {CSV} + {JSON}")
    return 0 if sidecar["all_checks_pass"] else 1


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