E-023-finite-height-assembly.py

{
  "experiment": "E-023",
  "task": "3.8 support: D-037 assembly (short-segment max, S_inf lower bound validation, margin/N0/N1)",
  "script": "research/numerics/scripts/E-023-finite-height-assembly.py",
  "checks": {
    "V1_sinf_lower_holds": true,
    "V2_margin_positive_convex": true,
    "V3_lambda_positive_3000": true
  },
  "key_values": {
    "M0_grid_max": 1.3323650421270834,
    "M0_lipschitz": 1.4075796709577815,
    "M0_use": 1.335180201468999,
    "C1": -1.1303307007539063,
    "C1_lower_form": -1.4189385332046727,
    "sinf_lower_min_slack": 9.520071987843544,
    "N0_measured_constants": 30217,
    "N0_traced_constants": 32871,
    "margin_N0_traced": 0.0795673057436943,
    "margin_at": {
      "32871": 0.0795673057436943,
      "100000": 287298.63095573196,
      "1000000": 5256686.503113639,
      "1000000000": 8938721534.500193,
      "1000000000000": 12396419838890.791,
      "1000000000000000000": 1.9304327085473796e+19,
      "8999999999994000000000001": 2.457962538414764e+26
    },
    "N1": 8999999999994000000000001,
    "N1_conservative_lagarias_window": 2.725978018024122e+21,
    "T0_platt_trudgian": 3000000000000.0
  },
  "elapsed_seconds": 56.2
}

"""E-023 — task 3.8 support: D-037 assembly computation.

Three components:
  S1. Short-segment max M0 = max |zeta'/zeta(-5/2+it)| for |t| <= 10.
      CLOSES THE D-036 GAP: the |R(s)| <= A log(t+2)+B0 bound was
      E-021-verified only for |t| >= 10; segment (III) passes within
      0.5 of the trivial pole at w = -2 (digamma pole), so the uniform
      bound fails near t = 0. Method: mpmath dps 30 grid (step 0.002)
      via the E-021 xi-route + sampled-Lipschitz correction x2 safety.
      (Only feeds an additive O(1) constant; slack is enormous.)
  S2. Validate the explicit S_inf lower bound used in D-037:
        S_inf(n) >= (1/2) n log n + (C1 - gamma/2) n - 4   (n >= 20)
      against E-021's EXACT S_inf (n <= 5000, validated 4.8e-25).
  S3. Assemble the full |S_f(n)| bound B(n) (D-036 constants, A = 0.9,
      B0 = 1.0, M0 from S1) and the margin
        margin(n) = Sinf_lower(n) - B(n) + 1;
      find N0 (margin > 0 and increasing); print margin at reference
      points and at N1 = (T0-1)^2, T0 = 3e12 (Platt-Trudgian 2021);
      conservative fallback N1* = T0^2/(4 log^2 T0) (Lagarias (1.16)).

Checks:
  V1. S_inf lower bound holds for all 20 <= n <= 5000 (exact values).
  V2. margin(n) > 0 for all n in [N0, 1e7] (grid, step 1) and margin
      convex (2nd differences > 0) on a log grid to N1 — with
      convexity + margin(N0) > 0 + margin'(N0) > 0 this gives
      margin(n) > 0 for ALL n in [N0, N1] (analytic in D-037).
  V3. lambda_n > 0 reconfirmed for n <= 3000 from E-021 exact CSV.

Outputs: research/numerics/data/E-023-assembly.csv (+ .json)
"""
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"
CSV = DATA / "E-023-assembly.csv"
JSON = DATA / "E-023-assembly.json"

GAMMA = 0.577215664901532860606512090082402431
C1 = (GAMMA - 1 - np.log(2 * np.pi)) / 2          # -1.1303...
C1L = C1 - GAMMA / 2                              # -1.4189... (lower-bound form)
A_USE, B0_USE = 0.9, 1.0
T0_PLATT = 3.0e12


# ---------------------------------------------------------------- S1
def zeta_logder(s):
    """zeta'/zeta via xi'/xi minus explicit terms (E-021 route, dps 30)."""
    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 short_segment_max():
    mp.mp.dps = 30
    ts = np.arange(-10.0, 10.0 + 1e-12, 0.002)
    vals = np.array([abs(zeta_logder(complex(-2.5, t))) for t in ts])
    gmax = float(vals.max())
    # sampled Lipschitz constant (max adjacent slope), x2 safety
    lips = float(np.max(np.abs(np.diff(vals))) / 0.002)
    m0 = gmax + 0.002 * lips          # (h/2)*L would be 0.001*L; x2
    return gmax, lips, m0


# ------------------------------------------------------- D-036 pieces
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, a, b0):
    return 3 * C0(T) ** 2 + a * np.log(T + 2) + b0


def B_bound(n, m0, a=A_USE, b0=B0_USE):
    """Full |S_f(n)| upper bound (D-036, variant-B contour, T = sqrt(n)+1)."""
    n = float(n)
    T = np.sqrt(n) + 1.0
    twoN = 2 * (F_count(T) + S_bound(T))
    I3_long = 4 * T * ((2 / 3) * C0(T) + a * np.log(T + 2) + b0)
    I3_short = 40.0 * m0                      # |k|<=2, length 20, max m0
    I24 = 2 * (np.e ** 4 + 1) * (7 * H_bound(T, a, b0) + 3.3)
    # I1 in log space: log I1 = log(2 sqrt(n) (e^{sqrt(n)/2}+1) 0.5699)
    #   - (2 sqrt(n) - 2) log 2  ->  hugely negative for sqrt(n) > 60
    r = np.sqrt(n)
    logI1 = np.log(2 * r * 0.5699) + r / 2 - (2 * r - 2) * np.log(2)
    I1 = np.exp(logI1) if logI1 > -700 else 0.0
    triv = 1.0
    return twoN + I3_long + I3_short + I24 + I1 + triv


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


def margin(n, m0, a=A_USE, b0=B0_USE):
    return sinf_lower(n) - B_bound(n, m0, a, b0) + 1.0


def find_N0(m0, a, b0):
    ngrid = np.arange(20, 400001)
    marg = np.array([margin(n, m0, a, b0) for n in ngrid])
    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
    return N0, cont, float(marg[ngrid == N0][0]) if N0 > 0 else float("nan")


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()

    print("S1: short-segment max ...", flush=True)
    gmax, lips, m0 = short_segment_max()
    print(f"  grid max = {gmax:.4f}, Lipschitz {lips:.1f}, "
          f"M0 = {m0:.4f}", flush=True)

    print("S2: S_inf lower-bound validation vs E-021 exact ...", flush=True)
    import csv as csvmod
    with open(E021) as f:
        rows = list(csvmod.DictReader(f))
    ns = np.array([int(r["n"]) for r in rows])
    sinf = np.array([float(r["S_inf"]) for r in rows])
    lam = np.array([float(r["lambda_n"]) for r in rows])
    mask = ns >= 20
    slack = sinf[mask] - sinf_lower(ns[mask])
    v1 = bool(np.all(slack > 0))
    print(f"  V1: min slack = {slack.min():.4f} at n = "
          f"{ns[mask][slack.argmin()]} (bound holds: {v1})", flush=True)
    v3 = bool(np.all(lam[ns <= 3000] > 0))
    print(f"  V3: lambda_n > 0 for n <= 3000 reconfirmed: {v3}", flush=True)

    print("S3: margin assembly ...", flush=True)
    # Two constant sets: measured (0.9, 1.0) and traced (2.0, 2.0, the
    # D-036 first-estimate analytic assembly).
    N0_meas, cont_meas, margN0_meas = find_N0(m0, 0.9, 1.0)
    N0_tr, cont_tr, margN0_tr = find_N0(m0, 2.0, 2.0)
    N0, cont = N0_tr, cont_tr
    print(f"  N0 (measured 0.9/1.0) = {N0_meas} (contiguous {cont_meas})")
    print(f"  N0 (traced   2.0/2.0) = {N0_tr} (contiguous {cont_tr})")

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

    N1 = int(3_000_000_000_000 - 1) ** 2
    N1_cons = T0_PLATT ** 2 / (4 * np.log(T0_PLATT) ** 2)
    refs = {int(n): float(margin(n, m0, 2.0, 2.0)) for n in
            [N0, 100000, 10 ** 6, 10 ** 9, 10 ** 12, 10 ** 18, N1]}
    print(f"  N0 = {N0} (traced; contiguous to 4e5: {cont}, "
          f"convex: {conv})")
    print(f"  margin(N0) = {refs[N0]:.1f}")
    print(f"  N1 = (T0-1)^2 = {N1:.6e};  N1* (Lagarias window) = "
          f"{N1_cons:.6e}")
    print(f"  margin(N1) = {refs[N1]:.3e}")

    v2 = bool(N0 > 0 and cont and conv and cont_meas)

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

    meta = {
        "experiment": "E-023",
        "task": "3.8 support: D-037 assembly (short-segment max, S_inf "
                "lower bound validation, margin/N0/N1)",
        "script": "research/numerics/scripts/E-023-finite-height-assembly.py",
        "checks": {"V1_sinf_lower_holds": v1,
                   "V2_margin_positive_convex": v2,
                   "V3_lambda_positive_3000": v3},
        "key_values": {
            "M0_grid_max": gmax, "M0_lipschitz": lips, "M0_use": m0,
            "C1": C1, "C1_lower_form": C1L,
            "sinf_lower_min_slack": float(slack.min()),
            "N0_measured_constants": N0_meas,
            "N0_traced_constants": N0_tr,
            "margin_N0_traced": margN0_tr,
            "margin_at": {str(k): float(v) for k, v in refs.items()},
            "N1": N1, "N1_conservative_lagarias_window": N1_cons,
            "T0_platt_trudgian": T0_PLATT,
        },
        "elapsed_seconds": round(time.time() - t0, 1),
    }
    JSON.write_text(json.dumps(meta, indent=2))
    print(f"  wrote {CSV.name}, {JSON.name}; "
          f"elapsed {meta['elapsed_seconds']}s")
    return 0 if (v1 and v2 and v3) else 1


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