E-034-dogleg-validation.py

{
  "experiment": "E-034",
  "title": "D-074 zero-free dog-leg repair DL-0: exact identity H_low = H_high - V_right + V_left verified; piece majorants with slack; the temporary dog-leg Bn_DL keeps margin > 0 to N1 = (3e12-1)^2",
  "session": "040",
  "thread": "ii-prime-6-DL",
  "design_doc": "D-074",
  "checks": {
    "V1_C_zeta3_certified": true,
    "V2_kernel_charge_certified": true,
    "V3_q_certified": true,
    "V4_identity_exact": true,
    "V5_piece_slack": true,
    "V6_fold_margin_positive": true,
    "VA_mpmath_vs_dirichlet": true
  },
  "all_checks_pass": true,
  "dps": 50,
  "constants": {
    "C_zeta3_partial": "0.164822681359",
    "C_zeta3_tail": "8.502e-9",
    "C_zeta3_charged": "0.5",
    "e^{5/2}+1": "13.182494",
    "q": "0.77981175",
    "val_max_disc": "4.97e-10",
    "foldfit_max": "0.3986",
    "margin_N1": "2.4579624e+26"
  },
  "notes": "Identity (*): H_low = H_high - V_right + V_left at T_cut = sqrt(n)+eps (zero-avoiding), U = 1.5*T_cut, X_R = 2*sqrt(n). DogLegRight(n) = VLeftBound + 5.94 + IIMaj/4 replaces the broken 3.3 (NOT multiplied by the mid kernel charge e^4+1)."
}

#!/usr/bin/env python3
"""E-034 — D-074 zero-free dog-leg repair, sub-stage DL-0 (directive §11):
numerical validation BEFORE any Lean work.

The dog-leg: keep T_cut = sqrt(n) + eps (zero-avoiding); replace the right
horizontal tail H_low = int_3^{2 sqrt(n)} Jlag(sigma + i T_cut) d sigma by
    H_low = H_high - V_right + V_left       (exact identity (*))
    H_high  = int_3^{2 sqrt(n)} Jlag(sigma + i U) d sigma,  U = 1.5 * T_cut
    V_right = i int_{T_cut}^{U} Jlag(2 sqrt(n) + i t) dt
    V_left  = i int_{T_cut}^{U} Jlag(3 + i t) dt
with Jlag n s = ((s/(s-1))^n - 1) * (-zeta'/zeta(s)) (the Lean convention,
LiSides.lean:55-58). The whole dog-leg rectangle is at Re s >= 3 > 1
(unconditionally zero-free).

CHECKS:
  V1: C_zeta3 certified: partial Dirichlet sum + certified tail brackets
      -zeta'/zeta(3) (mpmath cross-check); charged at the upper edge.
  V2: e^{5/2}+1 <= 13.19 certified (kernel charge on V_left).
  V3: q = e^{4/9}/2 <= 0.78 < 1 certified (the stepping condition).
  V4 (identity): |H_low - (H_high - V_right + V_left)| tiny at
      n = 1e3, 1e4, 1e5 (two trapezoid grids for convergence).
  V5 (slack): measured |V_left| <= VLeftBound, |H_high| <= 5.94,
      |V_right| <= VRightBound, |H_low| <= DogLegRight (ratios reported).
  V6 (fold): the temporary dog-leg Bn_DL(n) keeps margin > 0 on
      [8e4, N1], N1 = (3e12-1)^2 (full fold, base region, endpoint).

Method: dps 50; -zeta'/zeta via certified-tail Dirichlet series for
sigma >= 5 and mpmath for sigma in [3,5) (validated against it at
sigma = 3, 4, 5 — the E-032 architecture); vertical pieces by mpmath
(validated at sigma = 3); trapezoid h = 0.05 / 0.02 (horizontals),
h = 0.05 / 0.025 (verticals). Deterministic; cache via sidecar.
"""
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-034-dogleg-validation.csv"
JSON = DATA / "E-034-dogleg-validation.json"

SIEVE_MAX = 25000
TOL_DIR = mp.mpf("1e-11")
NS = [10 ** 3, 10 ** 4, 10 ** 5]
pi, e = mp.pi, mp.e

_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 tail_cert(sig, K):
    if K <= 3:
        return mp.inf
    u = mp.mpf(K - 1)
    return (mp.log(u) / (sig - 1) + 1 / (sig - 1) ** 2) * u ** (1 - sig)


def choose_K(sig, tol):
    K = 4
    while tail_cert(sig, K) >= tol:
        K += 1
        if K > SIEVE_MAX:
            raise RuntimeError(f"K exceeds sieve at sig={sig}")
    return K


def dser(sig, T, tol=TOL_DIR):
    K = choose_K(sig, tol)
    w = mp.mpc(sig, T)
    s = mp.mpc(0)
    for k in range(2, K + 1):
        if LAM[k]:
            s += LAM[k] * mp.exp(-w * LOGK[k])
    return s, K


def mp_neg_ld(sig, T):
    w = mp.mpc(sig, T)
    return -mp.zeta(w, derivative=1) / mp.zeta(w)


def neg_ld(sig, T):
    if sig < 5:
        return mp_neg_ld(sig, T), 0
    return dser(sig, T)


def jlag(n, sig, T):
    s = mp.mpc(sig, T)
    kf = (s / (s - 1)) ** n - 1
    d, _ = neg_ld(sig, T)
    return kf * d


def jlag_mp(n, sig, T):
    s = mp.mpc(sig, T)
    return ((s / (s - 1)) ** n - 1) * mp_neg_ld(sig, T)


def zeros_near(T0, half=3.0):
    j_approx = int((T0 / (2 * mp.pi)) * mp.log(T0 / (2 * mp.pi * mp.e)))
    out = []
    for j in range(max(1, j_approx - 30), j_approx + 30):
        g = mp.zetazero(j).imag
        if abs(g - T0) <= half:
            out.append(g)
    return out


def choose_eps(T0):
    zs = zeros_near(T0)
    best = (-1, 0.0)
    for i in range(20):
        eps = i * 0.05
        cl = min((abs(T0 + eps - g) for g in zs), default=mp.mpf(3))
        if cl > best[0]:
            best = (cl, eps)
    return best[1], best[0]


def H_piece(n, y, a, b, h):
    """int_a^b Jlag(sigma + i y) d sigma (mpmath < 5, Dirichlet >= 5)."""
    N = int(mp.ceil((b - a) / h))
    he = (b - a) / N
    tot = mp.mpc(0)
    for i in range(N + 1):
        sig = a + i * he
        f = jlag_mp(n, sig, y) if sig < 5 else jlag(n, sig, y)
        wt = mp.mpf("0.5") if i in (0, N) else mp.mpf(1)
        tot += wt * f
    return tot * he


def V_piece(n, x, a, b, h):
    """i int_a^b Jlag(x + i t) dt (mpmath throughout; validated at x = 3)."""
    N = int(mp.ceil((b - a) / h))
    he = (b - a) / N
    tot = mp.mpc(0)
    for i in range(N + 1):
        t = a + i * he
        f = jlag_mp(n, x, t)
        wt = mp.mpf("0.5") if i in (0, N) else mp.mpf(1)
        tot += wt * f
    return 1j * tot * he


# ------------------------------------------------- fold machinery (E-030 forms)
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 IIMaj(n):
    return 2 * (mp.sqrt(n) + 1) * (mp.e ** (mp.sqrt(n) / 2) + 1) * mp.mpf("0.57") * 2 ** (-(2 * mp.sqrt(n) - 2))


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


def NI(T, a):
    c_clear = 1 / (9 * C0hat(a + 2))
    return Z2(T) * 2 * mp.log(mp.mpf("16.5") / c_clear)


# candidate certified constants (V1-V3 below)
# C_ZETA3: the LEAN-CHARGED constant 0.5 (DL-2, dirichletC3 <= 1/2 via the
# termwise 3-vs-2 comparison + dirichletC2 <= 1). The tight analytic value
# -zeta'/zeta(3) = 0.16482 (V1 certifies it) and the fold is indifferent
# to the constant (margins stay > 1.3e4 at the tightest base point), so the
# cheap certified 0.5 is charged; the 0.165/0.17 upgrades are documented
# in D-074 sec.4 as the tightening path if ever needed.
C_ZETA3 = mp.mpf("0.5")
E52P1 = mp.mpf("13.19")          # e^{5/2}+1 <= 13.19
Q34 = mp.mpf("0.78")             # e^{4/9}/2 <= 0.78
# HHIGH: the LEAN-CHARGED high-horizontal bound 8.38 (DL-3,
# hIntegral_Jlag_dogleg_norm_le: the single-exponential-envelope form
# 4*e^{4/3}/(2c) + 1/(2 log 2) with c = log 2 - 4/9 >= 0.2487 certified,
# dirichletC2 <= 1). The tight R(1.5) = 5.94 would need a dirichletC2 <= 0.57
# build (documented tightening path; the fold is indifferent: margins change
# by ~5 at the base points against margins > 1.3e4).
HHIGH = mp.mpf("8.38")


def VLeftBound(n, cz3=C_ZETA3):
    return ((mp.sqrt(n) + 1) / 2) * E52P1 * cz3


def VRightBound(n):
    # LEAN-PROVABLE charge IIMaj/4 (DL-6-ii, vIntegral_Jlag_right_le_IIMaj_quarter:
    # the DIRECT dirichlet bound at re = 2 sqrt(n) — log2-leading + 3^{-(2 sqrt(n)-2)}
    # tail with dirichletC2 <= 1 — gives ||V_right|| <= IIMaj/4 EXACTLY, NO
    # dirichletC2 <= 0.57 needed; the fold is indifferent: IIMaj ~ 1e-105).
    return IIMaj(n) / 4


def DogLegRight(n, cz3=C_ZETA3):
    return VLeftBound(n, cz3) + HHIGH + VRightBound(n)


def Bn_DL(n, cz3=C_ZETA3):
    """Temporary dog-leg Bn (D-074 sec.6): the session-038 honest Bn with
    the broken right-tail 3.3 replaced by DogLegRight (NOT multiplied by
    the mid kernel charge)."""
    T = 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 * (e ** 4 + 1) * (7 * Hhat(T))
            + 2 * DogLegRight(n, cz3)
            + IIMaj(n) + 1)


def honest_DL(n, cz3=C_ZETA3):
    T = mp.sqrt(n) + 1
    a = mp.sqrt(n)
    mid = 2 * (e ** 4 + 1) * (NI(T, a) + mp.mpf("5.5") * (mp.mpf("10.5") * mp.log(T + 2) + 48))
    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") + mid + 2 * DogLegRight(n, cz3) + IIMaj(n) + 1)


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

    rows = []

    # ---- V1: C_zeta3 certified ----
    K3 = 25000
    s3 = mp.mpf(0)
    for k in range(2, K3 + 1):
        if LAM[k]:
            s3 += LAM[k] / (mp.mpf(k) ** 3)
    tail3 = (2 * mp.log(K3 - 1) + 1) / (4 * (K3 - 1) ** 2)
    z3mp = mp_neg_ld(3, 0).real
    v1 = s3 <= z3mp <= s3 + tail3 and s3 + tail3 <= C_ZETA3
    rows.append({"check": "V1_C_zeta3", "n": "", "quantity": "-zeta'/zeta(3)",
                 "measured": mp.nstr(z3mp, 12), "bound_lo": mp.nstr(s3, 12),
                 "bound_hi": mp.nstr(s3 + tail3, 12), "charged": mp.nstr(C_ZETA3, 6),
                 "pass": str(bool(v1))})

    # ---- V2/V3: kernel charge + q certified ----
    e52 = mp.e ** (mp.mpf(5) / 2)
    v2 = e52 + 1 <= E52P1
    rows.append({"check": "V2_kernel_charge_e^{5/2}+1", "n": "", "quantity": "e^{5/2}+1",
                 "measured": mp.nstr(e52 + 1, 8), "bound_lo": "",
                 "bound_hi": mp.nstr(E52P1, 6), "charged": mp.nstr(E52P1, 6),
                 "pass": str(bool(v2))})
    e49 = mp.e ** (mp.mpf(4) / 9)
    q = e49 / 2
    v3 = q <= Q34 and e49 < 2
    rows.append({"check": "V3_q_certified", "n": "", "quantity": "q = e^{4/9}/2",
                 "measured": mp.nstr(q, 8), "bound_lo": mp.nstr(e49, 8),
                 "bound_hi": mp.nstr(Q34, 6), "charged": mp.nstr(Q34, 6),
                 "pass": str(bool(v3))})

    # ---- validation: mpmath vs Dirichlet at sigma = 3, 4, 5 ----
    val_max = mp.mpf(0)
    for n in NS:
        T0 = mp.sqrt(n)
        eps, _cl = choose_eps(T0)
        T = T0 + eps
        for sig in [3.0, 4.0, 5.0]:
            dm, _k = dser(mp.mpf(sig), T, tol=mp.mpf("1e-8"))
            dd = mp_neg_ld(sig, T)
            if abs(dm - dd) > val_max:
                val_max = abs(dm - dd)

    # ---- V4/V5: the dog-leg pieces + identity + slack ----
    v4, v5 = True, True
    for n in NS:
        sq = mp.sqrt(n)
        T0 = sq
        eps, clearance = choose_eps(T0)
        T = T0 + eps
        U = mp.mpf("1.5") * T
        XR = 2 * sq
        Hlow1 = H_piece(n, T, 3, XR, 0.05)
        Hlow2 = H_piece(n, T, 3, XR, 0.02)
        Hhigh1 = H_piece(n, U, 3, XR, 0.05)
        Hhigh2 = H_piece(n, U, 3, XR, 0.02)
        Vr1 = V_piece(n, XR, T, U, 0.05)
        Vr2 = V_piece(n, XR, T, U, 0.025)
        Vl1 = V_piece(n, 3, T, U, 0.05)
        Vl2 = V_piece(n, 3, T, U, 0.025)
        # identity: H_low = H_high - V_right + V_left
        resid = abs(Hlow2 - (Hhigh2 - Vr2 + Vl2))
        conv = max(abs(Hlow1 - Hlow2), abs(Hhigh1 - Hhigh2), abs(Vr1 - Vr2), abs(Vl1 - Vl2))
        tol = mp.mpf("1e-6") * max(mp.mpf(1), abs(Hlow2)) + 4 * conv
        ok4 = resid <= tol
        v4 = v4 and ok4
        # slack vs majorants
        vlb, vrb, dlr = VLeftBound(n), VRightBound(n), DogLegRight(n)
        ratios = {
            "V_left": abs(Vl2) / vlb, "H_high": abs(Hhigh2) / HHIGH,
            "V_right": abs(Vr2) / vrb, "H_low": abs(Hlow2) / dlr,
        }
        ok5 = all(r <= 1 for r in ratios.values())
        v5 = v5 and ok5
        for name, val, bnd, ratio in [
                ("V_left", abs(Vl2), vlb, ratios["V_left"]),
                ("H_high", abs(Hhigh2), HHIGH, ratios["H_high"]),
                ("V_right", abs(Vr2), vrb, ratios["V_right"]),
                ("H_low", abs(Hlow2), dlr, ratios["H_low"])]:
            rows.append({"check": f"V5_slack_{name}", "n": str(n), "quantity": name,
                         "measured": mp.nstr(val, 6), "bound_lo": "",
                         "bound_hi": mp.nstr(bnd, 6), "charged": mp.nstr(ratio, 4),
                         "pass": str(bool(ratio <= 1))})
        rows.append({"check": "V4_identity", "n": str(n), "quantity": "H_low - (H_high - V_right + V_left)",
                     "measured": mp.nstr(resid, 3), "bound_lo": mp.nstr(conv, 3),
                     "bound_hi": mp.nstr(tol, 3),
                     "charged": f"T={mp.nstr(T, 8)} eps={eps} U={mp.nstr(U, 8)}",
                     "pass": str(bool(ok4))})
        print(f"n={n}: identity resid={mp.nstr(resid, 2)} (tol {mp.nstr(tol, 2)}) | "
              f"V_left {mp.nstr(abs(Vl2), 4)}/{mp.nstr(vlb, 4)} | "
              f"H_high {mp.nstr(abs(Hhigh2), 4)}/{mp.nstr(HHIGH, 4)} | "
              f"V_right {mp.nstr(abs(Vr2), 3)}/{mp.nstr(vrb, 3)} | "
              f"H_low {mp.nstr(abs(Hlow2), 4)}/{mp.nstr(dlr, 4)}", flush=True)

    # ---- V6: the temporary fold ----
    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)]

    N1 = (mp.mpf("3e12") - 1) ** 2
    grid = [mp.mpf("8e4"), mp.mpf("1e5"), mp.mpf("2e5"), mp.mpf("5e5")] + \
        loggrid(mp.mpf("1e6"), N1, 25)
    v6 = True
    minmargin = None
    foldfit_max = mp.mpf(0)
    for nn in grid:
        T = mp.sqrt(nn) + 1
        a = 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 > foldfit_max:
            foldfit_max = fit
        marg = sInftyLower(nn) - Bn_DL(nn) + 1
        if minmargin is None or marg < minmargin:
            minmargin = marg
        v6 = v6 and (marg > 0) and (honest_DL(nn) <= Bn_DL(nn))
    for nn in [mp.mpf("8e4"), mp.mpf("1e5"), mp.mpf("1e8"), mp.mpf("1e12"), N1]:
        marg = sInftyLower(nn) - Bn_DL(nn) + 1
        rows.append({"check": "V6_margin_value", "n": mp.nstr(nn, 8), "quantity": "margin_DL",
                     "measured": mp.nstr(marg, 8), "bound_lo": "0",
                     "bound_hi": "", "charged": mp.nstr(Bn_DL(nn), 8),
                     "pass": str(bool(marg > 0))})
    rows.append({"check": "V6_margin_min_on_grid", "n": "grid", "quantity": "min margin",
                 "measured": mp.nstr(minmargin, 6), "bound_lo": "0", "bound_hi": "",
                 "charged": mp.nstr(foldfit_max, 4), "pass": str(bool(minmargin > 0 and foldfit_max < 1))})
    v6 = v6 and foldfit_max < 1

    allpass = v1 and v2 and v3 and v4 and v5 and v6 and val_max <= mp.mpf("1e-6")

    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-034",
        "title": "D-074 zero-free dog-leg repair DL-0: exact identity "
                 "H_low = H_high - V_right + V_left verified; piece majorants "
                 "with slack; the temporary dog-leg Bn_DL keeps margin > 0 to "
                 "N1 = (3e12-1)^2",
        "session": "040",
        "thread": "ii-prime-6-DL",
        "design_doc": "D-074",
        "checks": {"V1_C_zeta3_certified": bool(v1),
                   "V2_kernel_charge_certified": bool(v2),
                   "V3_q_certified": bool(v3),
                   "V4_identity_exact": bool(v4),
                   "V5_piece_slack": bool(v5),
                   "V6_fold_margin_positive": bool(v6),
                   "VA_mpmath_vs_dirichlet": bool(val_max <= mp.mpf("1e-6"))},
        "all_checks_pass": bool(allpass),
        "dps": 50,
        "constants": {"C_zeta3_partial": mp.nstr(s3, 12),
                      "C_zeta3_tail": mp.nstr(tail3, 4),
                      "C_zeta3_charged": mp.nstr(C_ZETA3, 6),
                      "e^{5/2}+1": mp.nstr(e52 + 1, 8),
                      "q": mp.nstr(q, 8),
                      "val_max_disc": mp.nstr(val_max, 3),
                      "foldfit_max": mp.nstr(foldfit_max, 4),
                      "margin_N1": mp.nstr(sInftyLower(N1) - Bn_DL(N1) + 1, 8)},
        "notes": "Identity (*): H_low = H_high - V_right + V_left at "
                 "T_cut = sqrt(n)+eps (zero-avoiding), U = 1.5*T_cut, X_R = 2*sqrt(n). "
                 "DogLegRight(n) = VLeftBound + 5.94 + IIMaj/4 replaces the broken 3.3 "
                 "(NOT multiplied by the mid kernel charge e^4+1).",
    }
    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)} VA {mp.nstr(val_max, 2)}")
    print(f"wrote {CSV} + {JSON}")
    return 0 if sidecar["all_checks_pass"] else 1


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