E-032-horizontal-true-integral.py

{
  "experiment": "E-032",
  "title": "F-D side thread 6-R-0: the TRUE complex right-tail horizontal integral of Lagarias Thm 6.1 at T = sqrt(n) vs the E-031 envelope",
  "session": "039",
  "thread": "ii-prime-6-R",
  "design_doc": "D-073",
  "finding_doc": "D-072",
  "checks": {
    "V0_envelope_reproduces_E031": true,
    "VA_mpmath_vs_dirichlet_le_1e-6": true,
    "VC_trapezoid_convergence": true,
    "V3_repair_envelope_collapses": true
  },
  "all_checks_pass": true,
  "dps": 50,
  "method": "-zeta'/zeta by truncated Dirichlet series (certified tail 1e-11) for sigma >= 5; mpmath validated against it at {3,3.5,4,4.5,5} for sigma in [3,5); trapezoid h=0.05/0.02.",
  "verdict_inputs": {
    "1000": {
      "Re_I_right": "2.1287636",
      "Im_I_right": "-6.8620183",
      "abs_I_right": "7.18463",
      "envelope_right": "50.8152",
      "abs_over_env": "0.1414",
      "abs_over_sqrtn": "0.2272",
      "abs_I_mid": "4.38071",
      "envelope_mid": "4.58947",
      "Re_I_full": "-2.2458021",
      "abs_I_full": "7.44101",
      "abs_I_right_c1.5": "1.12202",
      "abs_I_right_c2": "0.623246"
    },
    "10000": {
      "Re_I_right": "-11.887447",
      "Im_I_right": "-11.73172",
      "abs_I_right": "16.7016",
      "envelope_right": "10258.0",
      "abs_over_env": "0.001628",
      "abs_over_sqrtn": "0.167",
      "abs_I_mid": "5.22324",
      "envelope_mid": "7.41556",
      "Re_I_full": "-16.646154",
      "abs_I_full": "19.2052",
      "abs_I_right_c1.5": "1.14022",
      "abs_I_right_c2": "0.312063"
    },
    "100000": {
      "Re_I_right": "-13.538954",
      "Im_I_right": "22.464752",
      "abs_I_right": "26.2292",
      "envelope_right": "4.70658e+10",
      "abs_over_env": "5.573e-10",
      "abs_over_sqrtn": "0.08294",
      "abs_I_mid": "13.6968",
      "envelope_mid": "14.739",
      "Re_I_full": "-25.972869",
      "abs_I_full": "30.8894",
      "abs_I_right_c1.5": "0.999556",
      "abs_I_right_c2": "0.354185"
    }
  },
  "notes": "Verdict rule (D-073): |I_right| ~ envelope => per-side claim FALSE at T = sqrt(n); |I_right| << envelope => hidden cancellation (STOP and report)."
}

#!/usr/bin/env python3
"""E-032 — F-D side thread ii′-6-R, sub-stage 6-R-0 (deliverable D1):
the TRUE complex right-tail horizontal integral of Lagarias's Thm 6.1
contour at T = sqrt(n) (zero-avoiding) vs the E-031 envelope.

QUESTION (D-073 §0): is the per-side claim |I_II| = O(sqrt(n)) FALSE at
T = sqrt(n) — i.e. is the true complex integral
    I_right = int_{3}^{2 sqrt(n)} k_n(w-1) . (-zeta'/zeta(w)) d sigma,
    w = sigma + iT,
of the SAME order as the envelope int |...| d sigma (~5.9e10 at n=1e5,
E-031), or is it MUCH SMALLER via oscillation (hidden cancellation)?

METHOD (dps 50):
  * -zeta'/zeta(w) for sigma >= 5: TRUNCATED DIRICHLET SERIES
    sum_{k<=K} Lambda(k)/k^w with K chosen per sigma so the certified
    tail bound [log(K-1)/(sig-1) + 1/(sig-1)^2] . (K-1)^{1-sig} < 1e-11.
    Exact (this IS -zeta'/zeta for re > 1), no library zeta involved.
  * sigma in [3, 5): mpmath -zeta'/zeta (moderate args, trustworthy),
    VALIDATED against the Dirichlet series at sigma in {3, 3.5, 4, 4.5, 5}
    (tol 1e-8, K ~ 2.3e4 at sigma=3).
  * kernel k_n(w-1) = (1 + 1/(w-1))^n - 1, direct mpmath power.
  * T = c.sqrt(n) + eps, eps from a deterministic zero-avoiding scan
    (max-clearance grid point, zeros via mp.zetazero).
  * Composite trapezoid, two grids h = 0.05 / 0.02; the h-delta is the
    quadrature-error estimate. Integrand analytic, |phase'| <= ~2
    rad/unit, endpoints tiny vs peak => near-spectral convergence.

CHECKS:
  V0: the refined envelope at c=1 reproduces E-031 (n=1e4: ~1.05e4;
      n=1e5: ~5.9e10).
  VA: mpmath-vs-Dirichlet validation max discrepancy <= 1e-6.
  VC: trapezoid h-delta <= 1e-4 . envelope.
  V3: the repair collapses the envelope: Env(c=2)/Env(c=1) < 1e-6 (n=1e5).

VERDICT FIELDS (reported, not auto-decided): Re/Im/|I_right|, envelope,
|I|/envelope, and for c=1 also the mid part [-5/2, 3] (mpmath zeta) and
the full I_II = I_mid + I_right.
"""
import json
import math
import sys
from pathlib import Path

import mpmath as mp

mp.mp.dps = 50

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

SIEVE_MAX = 25000
TOL_DIR = mp.mpf("1e-11")     # Dirichlet tail certificate target (sigma >= 5)
TOL_SPOT = mp.mpf("1e-8")     # validation spot checks (sigma in {3..5})
NS = [10 ** 3, 10 ** 4, 10 ** 5]
CS = [1.0, 1.5, 2.0]

# ---------------------------------------------------------------- von Mangoldt
_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)   # Lambda(k) as mpf log p, else 0
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):
    """Certified bound on sum_{k>K} Lambda(k)/k^sig (sig >= 3, K >= 3):
    <= integral_{K-1}^inf log x / x^sig dx (summand decreasing)."""
    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):
    """-zeta'/zeta(sig + iT) via the truncated Dirichlet series, re > 1."""
    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):
    """Primary -zeta'/zeta: mpmath on [3,5), Dirichlet series on [5,inf)."""
    if sig < 5:
        return mp_neg_ld(sig, T), 0
    return dser(sig, T)


def kernel(n, sig, T):
    return (1 + 1 / mp.mpc(sig - 1, T)) ** n - 1


def integrand(n, sig, T):
    d, _ = neg_ld(sig, T)
    return kernel(n, sig, T) * d


def zeros_near(T0, half=3.0):
    """Zero ordinates within [T0-half, T0+half] (deterministic)."""
    j_approx = int((T0 / (2 * math.pi)) * math.log(T0 / (2 * math.pi * math.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):
    """Deterministic zero-avoiding eps in {0, .05, ..., .95}."""
    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 integrate(n, T, a, b, h, use_mid=False):
    """Trapezoid for I = int_a^b f d sigma and the envelope int |f|."""
    N = int(math.ceil((b - a) / h))
    he = (b - a) / N
    tot = mp.mpc(0)
    env = mp.mpf(0)
    for i in range(N + 1):
        sig = a + i * he
        if use_mid or sig < 5:
            f = kernel(n, sig, T) * mp_neg_ld(sig, T)
        else:
            f = integrand(n, sig, T)
        wt = mp.mpf("0.5") if i in (0, N) else mp.mpf(1)
        tot += wt * f
        env += wt * abs(f)
    return tot * he, env * he


def right_part(n, T, h):
    return integrate(n, T, 3.0, 2 * mp.sqrt(n), h)


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

    rows = []
    results = {}
    max_val_disc = mp.mpf(0)

    for n in NS:
        sq = mp.sqrt(n)
        for c in CS:
            T0 = c * sq
            eps, clearance = choose_eps(T0)
            T = T0 + eps
            # ---- validation (mpmath vs Dirichlet spot checks) ----
            val = {}
            for sig in [3.0, 3.5, 4.0, 4.5, 5.0]:
                dm, _k = dser(mp.mpf(sig), T, tol=TOL_SPOT)
                dd = mp_neg_ld(sig, T)
                disc = abs(dm - dd)
                val[sig] = disc
                if disc > max_val_disc:
                    max_val_disc = disc
            # ---- right part, two grids ----
            I1, E1 = right_part(n, T, 0.05)
            I2, E2 = right_part(n, T, 0.02)
            rec = {
                "n": n, "c": c, "T": T, "eps": eps, "clearance": clearance,
                "I_right": I2, "env": E2, "hdelta": abs(I1 - I2),
                "env_hdelta": abs(E1 - E2),
                "val_max": max(val.values()),
            }
            # ---- mid part + full side (c = 1 only) ----
            if c == 1.0:
                M1, ME1 = integrate(n, T, -2.5, 3.0, 0.05, use_mid=True)
                M2, ME2 = integrate(n, T, -2.5, 3.0, 0.02, use_mid=True)
                rec.update({"I_mid": M2, "env_mid": ME2,
                            "mid_hdelta": abs(M1 - M2),
                            "I_full": M2 + I2})
            results[(n, c)] = rec
            print(f"n={n} c={c}: T={mp.nstr(T,8)} eps={eps} "
                  f"|I_right|={mp.nstr(abs(I2),4)} env={mp.nstr(E2,4)} "
                  f"ratio={mp.nstr(abs(I2)/E2,3)} hdelta={mp.nstr(rec['hdelta'],2)}",
                  flush=True)
            rows.append({
                "check": "right_part", "n": str(n), "c": mp.nstr(c, 3),
                "T": mp.nstr(T, 10), "eps": str(eps),
                "clearance": mp.nstr(clearance, 4),
                "Re_I_right": mp.nstr(I2.real, 8),
                "Im_I_right": mp.nstr(I2.imag, 8),
                "abs_I_right": mp.nstr(abs(I2), 6),
                "envelope": mp.nstr(E2, 6),
                "abs_over_env": mp.nstr(abs(I2) / E2, 4),
                "Re_over_env": mp.nstr(abs(I2.real) / E2, 4),
                "abs_over_sqrtn": mp.nstr(abs(I2) / sq, 4),
                "hdelta": mp.nstr(rec["hdelta"], 3),
                "val_max_disc": mp.nstr(rec["val_max"], 3),
                "abs_I_mid": mp.nstr(abs(rec["I_mid"]), 6) if c == 1.0 else "",
                "envelope_mid": mp.nstr(rec["env_mid"], 6) if c == 1.0 else "",
                "Re_I_full": mp.nstr(rec["I_full"].real, 8) if c == 1.0 else "",
                "abs_I_full": mp.nstr(abs(rec["I_full"]), 6) if c == 1.0 else "",
            })

    # ---------------------------------------------------------------- checks
    r = results
    v0 = (mp.mpf("0.9e4") < r[(10 ** 4, 1.0)]["env"] < mp.mpf("1.2e4")) and \
         (mp.mpf("4e10") < r[(10 ** 5, 1.0)]["env"] < mp.mpf("8e10"))
    va = max_val_disc <= mp.mpf("1e-6")
    vc = all(rec["hdelta"] <= mp.mpf("1e-4") * rec["env"] for rec in r.values())
    v3 = r[(10 ** 5, 2.0)]["env"] / r[(10 ** 5, 1.0)]["env"] < mp.mpf("1e-6")

    DATA.mkdir(exist_ok=True)
    with CSV.open("w") as f:
        f.write(",".join(rows[0].keys()) + "\n")
        for row in rows:
            f.write(",".join(row.values()) + "\n")

    verdict = {}
    for n in NS:
        rec = r[(n, 1.0)]
        verdict[str(n)] = {
            "Re_I_right": mp.nstr(rec["I_right"].real, 8),
            "Im_I_right": mp.nstr(rec["I_right"].imag, 8),
            "abs_I_right": mp.nstr(abs(rec["I_right"]), 6),
            "envelope_right": mp.nstr(rec["env"], 6),
            "abs_over_env": mp.nstr(abs(rec["I_right"]) / rec["env"], 4),
            "abs_over_sqrtn": mp.nstr(abs(rec["I_right"]) / mp.sqrt(n), 4),
            "abs_I_mid": mp.nstr(abs(rec["I_mid"]), 6),
            "envelope_mid": mp.nstr(rec["env_mid"], 6),
            "Re_I_full": mp.nstr(rec["I_full"].real, 8),
            "abs_I_full": mp.nstr(abs(rec["I_full"]), 6),
            "abs_I_right_c1.5": mp.nstr(abs(r[(n, 1.5)]["I_right"]), 6),
            "abs_I_right_c2": mp.nstr(abs(r[(n, 2.0)]["I_right"]), 6),
        }
    sidecar = {
        "experiment": "E-032",
        "title": "F-D side thread 6-R-0: the TRUE complex right-tail horizontal "
                 "integral of Lagarias Thm 6.1 at T = sqrt(n) vs the E-031 envelope",
        "session": "039",
        "thread": "ii-prime-6-R",
        "design_doc": "D-073",
        "finding_doc": "D-072",
        "checks": {"V0_envelope_reproduces_E031": bool(v0),
                   "VA_mpmath_vs_dirichlet_le_1e-6": bool(va),
                   "VC_trapezoid_convergence": bool(vc),
                   "V3_repair_envelope_collapses": bool(v3)},
        "all_checks_pass": bool(v0 and va and vc and v3),
        "dps": 50,
        "method": "-zeta'/zeta by truncated Dirichlet series (certified tail "
                  "1e-11) for sigma >= 5; mpmath validated against it at "
                  "{3,3.5,4,4.5,5} for sigma in [3,5); trapezoid h=0.05/0.02.",
        "verdict_inputs": verdict,
        "notes": "Verdict rule (D-073): |I_right| ~ envelope => per-side claim "
                 "FALSE at T = sqrt(n); |I_right| << envelope => hidden "
                 "cancellation (STOP and report).",
    }
    with JSON.open("w") as f:
        json.dump(sidecar, f, indent=2)
    print(f"V0 envelope~E031: {bool(v0)} | VA validation: {bool(va)} "
          f"(max {mp.nstr(max_val_disc,2)}) | VC convergence: {bool(vc)} | "
          f"V3 repair: {bool(v3)}")
    print(f"wrote {CSV} + {JSON}")
    return 0 if sidecar["all_checks_pass"] else 1


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