E-019-zeroflow-a12.py

{
  "experiment": "E-019",
  "task": "2.8(ii): backend v2 validation + zero-flow at a=12",
  "script": "research/numerics/scripts/E-019-zeroflow-a12.py",
  "parameters": {
    "a_val": 10.0,
    "a_new": 12.0,
    "h": 0.01,
    "lambda": -1.0,
    "z_max": 300.0
  },
  "checks": {
    "VA1_backend_matches": true,
    "VA2_pipeline_a10": true,
    "B3_no_new_complex_zeros": true
  },
  "report": {
    "validation": {
      "VA1_max_dev": 2.2246968001127243e-07,
      "VA2_Q_dev": 4.3871651103644105e-08,
      "VA2_z1_dev": 6.942222297023903e-09,
      "VA2_count_diff": 0
    },
    "B1_zeroflow": {
      "a10": {
        "n": 954,
        "z1": 0.44959701245436284,
        "z1_times_a": 4.4959701245436285,
        "count_vs_lattice": 954.9296585513721,
        "gap_mean": 0.31421194329970964,
        "gap_cv": 0.19770812588191275
      },
      "a12": {
        "n": 1145,
        "z1": 0.3746087478472408,
        "z1_times_a": 4.495304974166889,
        "count_vs_lattice": 1145.9155902616465,
        "gap_mean": 0.26188092266824087,
        "gap_cv": 0.18590651886845952
      }
    },
    "B2_gamma_locking": {
      "a10": {
        "min": 0.00020978025152373903,
        "zero_idx": 245,
        "gamma_idx": 20,
        "zero": 77.14504984912632,
        "gamma": 77.1448400688748,
        "median": 48.899737666799254
      },
      "a12": {
        "min": 5.2439703061679666e-05,
        "zero_idx": 95,
        "gamma_idx": 3,
        "zero": 25.010805140442628,
        "gamma": 25.01085758014569,
        "median": 48.79267672793192
      }
    },
    "B3_strip_winding": {
      "a10": {
        "1.0": {
          "winding": 253.99999999999946,
          "real_inside": 252,
          "excess": 1.99999999999946
        },
        "0.1": {
          "winding": 253.99999999999946,
          "real_inside": 252,
          "excess": 1.99999999999946
        }
      },
      "a12": {
        "1.0": {
          "winding": 305.9999999999994,
          "real_inside": 304,
          "excess": 1.9999999999993747
        },
        "0.1": {
          "winding": 305.9999999999993,
          "real_inside": 304,
          "excess": 1.9999999999993179
        }
      }
    },
    "B4_envelope": {
      "median_M_over_R0_y0.3": 0.006526107920803935
    }
  },
  "elapsed_seconds": 672.2,
  "timestamp_utc": "2026-08-18T19:32:51.149654+00:00"
}

"""E-019 — task 2.8(ii): backend v2 validation (A) + zero-flow at a=12 (B).

Sub-stage A (validation gate, all must pass before a=12 is trusted):
  V-A1. prime_backend_v2 vs E-016 in-memory cumsums: max |dev| on 1e4
        random t in (0,20] < 1e-6 abs (float64 summation-reordering
        noise floor ~2e-7 measured; both backends carry the same
        rounding level, so agreement at this level = exact).
  V-A2. full pipeline at a=10 with the v2 backend: stiffness matrix Q
        vs E16.build(10) entrywise max dev < 1e-8; located zeros at
        a=10 match the E-016 family (z1 dev < 1e-6; count equal).

Sub-stage B (measurements at a=12, theta=pi):
  B1. zero-flow: located real zeros (0,300]: count vs aT/pi, z1*a law,
      gap statistics vs a=10 (inflow continues? any bounce?).
  B2. gamma-locking metric: min |z_j - gamma_k| over first 30 Riemann
      zeros (mpmath zetazero, safe per STATE) vs a=10 value.
  B3. strip winding counts: W-zeros in |Re z| <= 40, |Im z| <= y0 for
      y0 in {1, 0.1} (argument principle; detects any non-real
      deviation beyond the known origin pair).
  B4. envelope spot check: median |M|/|R_0| at y0 = 0.3 (E-018 family).

Outputs: research/numerics/data/E-019-zeroflow-a12.csv (+ .json)
"""
import importlib.util
import json
import sys
import time
from datetime import datetime, timezone
from pathlib import Path

import mpmath as mp
import numpy as np
import scipy.optimize

REPO = Path(__file__).resolve().parents[3]
DATA = REPO / "research" / "numerics" / "data"
CSV = DATA / "E-019-zeroflow-a12.csv"
JSON = DATA / "E-019-zeroflow-a12.json"

_spec = importlib.util.spec_from_file_location(
    "E16", Path(__file__).with_name("E-016-herglotz-fine-structure.py"))
E16 = importlib.util.module_from_spec(_spec)
_spec.loader.exec_module(E16)                     # runs the e^20 sieve

sys.path.insert(0, str(Path(__file__).parent))
from prime_backend_v2 import PrimeBackendV2       # noqa: E402

H = E16.H_MAIN
LAM = E16.LAM_SHIFT
A_VAL = 10.0
A_NEW = 12.0
Z_MAX = 300.0
N_GRID = 40000
mp.mp.dps = 40

XN, XW = E16.XN, E16.XW


def collect_t_points(a, h):
    """Replicates E-016 S_full/quad segment structure (no g evaluation)
    to enumerate every quadrature abscissa t in (0, 2a)."""
    n_elem = int(round(2 * a / h))
    pts = []

    def add_quad(lo, hi, peak):
        n_lo = int(np.ceil(np.exp(lo)))
        n_hi = int(np.floor(np.exp(hi)))
        if n_hi - n_lo > 64:
            for aa, bb in ((lo, peak), (peak, hi)):
                if bb - aa > 0:
                    mid, half = 0.5 * (aa + bb), 0.5 * (bb - aa)
                    pts.extend(mid + half * XN)
            return
        segs = sorted({lo, peak, hi} |
                      {float(np.log(n))
                       for n in range(max(n_lo, 2), n_hi + 1)
                       if lo < np.log(n) < hi})
        for p0, p1 in zip(segs, segs[1:]):
            if p1 - p0 > 0:
                mid, half = 0.5 * (p0 + p1), 0.5 * (p1 - p0)
                pts.extend(mid + half * XN)

    for k in range(n_elem):
        if k == 0:
            add_quad(0.0, h, 0.0)
        else:
            add_quad((k - 1) * h, (k + 1) * h, k * h)
    return np.array(pts)


def locate_zeros(cp, x, h, z_max=Z_MAX, n_grid=N_GRID):
    zg = np.linspace(1e-3, z_max, n_grid)
    Vpr, _ = E16.V_and_dV(zg, cp, x, h)
    g = zg * Vpr.imag - Vpr.real

    def g_direct(zz):
        V = E16.hat_FT(np.array([zz]), h)[0] * complex(
            np.exp(1j * zz * x) @ cp)
        return zz * V.imag - V.real

    idx = np.flatnonzero(np.sign(g[:-1]) * np.sign(g[1:]) < 0)
    zeros = []
    for i in idx:
        try:
            zeros.append(scipy.optimize.brentq(g_direct, zg[i], zg[i + 1],
                                               xtol=1e-13, rtol=1e-13))
        except ValueError:
            pass
    return np.array(zeros)


def winding_rect(cp, cm, x, h, half_re, half_im, n_per_edge=4000):
    """Zeros of W(theta=pi) inside |Re z|<=half_re, |Im z|<=half_im via
    the argument principle (direct evaluation, e^{+-a*y0} range OK at
    y0 <= 1, a <= 12)."""
    edges = []
    xr = np.linspace(-half_re, half_re, n_per_edge)
    yr = np.linspace(-half_im, half_im, max(8, n_per_edge // 10))
    edges.append(xr - 1j * half_im)          # bottom, left->right
    edges.append(half_re + 1j * yr)          # right, bottom->top
    edges.append(xr[::-1] + 1j * half_im)    # top, right->left
    edges.append(-half_re + 1j * yr[::-1])   # left, top->bottom
    total = 0.0
    for z in edges:
        Vp, _ = E16.V_and_dV(z, cp, x, h)
        Vm, _ = E16.V_and_dV(z, cm, x, h)
        W = (z - 1j) * Vp - (z + 1j) * Vm
        total += np.diff(np.unwrap(np.angle(W))).sum()
    return total / (2 * np.pi)


def solve_ab(a, h):
    x, hh, Q, M, b = E16.build(a, h)
    n = len(x)
    cp, T = E16.solve_v(Q, M, b[+1.0], LAM, n)
    cm, _ = E16.solve_v(Q, M, b[-1.0], LAM, n)
    np_ = float(cp @ T @ cp)
    nm = float(cm @ T @ cm)
    return x, hh, cp, cm * np.sqrt(np_ / nm), Q


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

    # ---------------- backend v2: collect queries & build -------------
    print("collecting quadrature t-points (a=10 & a=12) ...", flush=True)
    rng = np.random.default_rng(19)
    t_rand = np.sort(rng.uniform(0.01, 20.0, 10000))
    tq = np.concatenate([collect_t_points(A_VAL, H),
                         collect_t_points(A_NEW, H), t_rand])
    print(f"  {tq.size} points; building backend to e^24 ...", flush=True)
    CACHE = DATA / "E-019-backend-cache.npz"
    be = PrimeBackendV2(tq)
    if CACHE.exists() and "--force" not in sys.argv:
        be.load_cache(CACHE)
        print("  backend loaded from cache", flush=True)
    else:
        be.build(verbose=True)
        be.save_cache(CACHE)

    # ---------------- A: validation -----------------------------------
    print("V-A1: backend vs E-016 in-memory cumsums ...", flush=True)
    p_ref = E16.prime_sum_vec(t_rand)
    p_new = be.prime_sum(t_rand)
    va1 = float(np.max(np.abs(p_ref - p_new)))
    print(f"  V-A1 max |dev|: {va1:.2e}", flush=True)

    print("V-A2: full pipeline at a=10 (v2 backend) ...", flush=True)
    orig_psv = E16.prime_sum_vec
    # reference build with the ORIGINAL backend first
    x10r, hh10r, cp10r, cm10r, Q_ref = solve_ab(A_VAL, H)
    zeros10r = locate_zeros(cp10r, x10r, hh10r)
    # rebuild with v2 backend (cache must be cleared — keyed on (a,h))
    E16._BUILD.clear()
    E16.prime_sum_vec = be.prime_sum
    x10, hh10, cp10, cm10, Q_new = solve_ab(A_VAL, H)
    E16.prime_sum_vec = orig_psv
    va2_q = float(np.max(np.abs(Q_new - Q_ref)))
    zeros10 = locate_zeros(cp10, x10, hh10)
    z1_10 = float(zeros10[0])
    va2_z1 = abs(z1_10 - float(zeros10r[0]))
    va2_cnt = abs(len(zeros10) - len(zeros10r))
    print(f"  V-A2 Q dev {va2_q:.2e}; z1 dev {va2_z1:.2e}; "
          f"count diff {va2_cnt}", flush=True)
    report["validation"] = {"VA1_max_dev": va1, "VA2_Q_dev": va2_q,
                            "VA2_z1_dev": va2_z1, "VA2_count_diff": va2_cnt}

    # ---------------- B: a=12 measurements ----------------------------
    print("B: a=12 build + zeros + winding ...", flush=True)
    E16.prime_sum_vec = be.prime_sum
    x12, hh12, cp12, cm12, _ = solve_ab(A_NEW, H)
    E16.prime_sum_vec = orig_psv
    zeros12 = locate_zeros(cp12, x12, hh12)

    # B1 zero-flow
    gaps10 = np.diff(zeros10)
    gaps12 = np.diff(zeros12)
    b1 = {
        "a10": {"n": int(len(zeros10)), "z1": z1_10,
                "z1_times_a": z1_10 * A_VAL,
                "count_vs_lattice": A_VAL * Z_MAX / np.pi,
                "gap_mean": float(gaps10.mean()),
                "gap_cv": float(gaps10.std() / gaps10.mean())},
        "a12": {"n": int(len(zeros12)), "z1": float(zeros12[0]),
                "z1_times_a": float(zeros12[0]) * A_NEW,
                "count_vs_lattice": A_NEW * Z_MAX / np.pi,
                "gap_mean": float(gaps12.mean()),
                "gap_cv": float(gaps12.std() / gaps12.mean())},
    }
    print(f"  B1 a=12: n={b1['a12']['n']} (lattice "
          f"{b1['a12']['count_vs_lattice']:.0f}) z1*a="
          f"{b1['a12']['z1_times_a']:.4f} gap_cv={b1['a12']['gap_cv']:.4f}",
          flush=True)

    # B2 gamma-locking
    gammas = np.array([float(mp.zetazero(k).imag) for k in range(1, 31)])
    def lock_metric(zeros):
        d = np.abs(zeros[:, None] - gammas[None, :])
        j, k = np.unravel_index(np.argmin(d), d.shape)
        dmin = np.min(d, axis=1)
        return (float(d[j, k]), int(j + 1), int(k + 1),
                float(zeros[j]), float(gammas[k]),
                float(np.median(dmin)))
    lm10 = lock_metric(zeros10)
    lm12 = lock_metric(zeros12)
    b2 = {"a10": {"min": lm10[0], "zero_idx": lm10[1], "gamma_idx": lm10[2],
                  "zero": lm10[3], "gamma": lm10[4], "median": lm10[5]},
          "a12": {"min": lm12[0], "zero_idx": lm12[1], "gamma_idx": lm12[2],
                  "zero": lm12[3], "gamma": lm12[4], "median": lm12[5]}}
    print(f"  B2 lock: a=10 min {lm10[0]:.5f} (z_{lm10[1]}={lm10[3]:.4f} "
          f"vs g_{lm10[2]}={lm10[4]:.4f}); a=12 min {lm12[0]:.5f} "
          f"(z_{lm12[1]}={lm12[3]:.4f} vs g_{lm12[2]}={lm12[4]:.4f})",
          flush=True)

    # B3 strip winding (a=12 and a=10 for comparison)
    b3 = {}
    for tag, cp, cm, x, hh in (("a10", cp10, cm10, x10, hh10),
                               ("a12", cp12, cm12, x12, hh12)):
        real_in = {}
        for y0 in (1.0, 0.1):
            w = winding_rect(cp, cm, x, hh, 40.0, y0)
            zz = zeros10 if tag == "a10" else zeros12
            n_real = int(2 * np.sum(zz <= 40.0))
            real_in[str(y0)] = {"winding": w, "real_inside": n_real,
                                "excess": w - n_real}
        b3[tag] = real_in
    print(f"  B3: {b3}", flush=True)

    # B4 envelope spot check at y0 = 0.3 (E-018 family), a=12
    def R0_exact(a, th, z):
        z = np.asarray(z, dtype=complex)
        t1 = -2j * (np.sinh(a * (1 + 1j * z))
                    + np.exp(1j * th) * np.sinh(a * (1j * z - 1)))
        t2 = (2j * np.sin(a * z) / z) * ((z - 1j) - np.exp(1j * th)
                                         * (z + 1j))
        return t1 + t2
    xg = np.linspace(-20, 20, 41) + 0.3j
    xm = x12.min()
    u = x12 - xm
    E = np.exp(1j * np.outer(xg, u))
    hat = E16.hat_FT(xg, hh12)
    W = np.exp(1j * xm * xg) * ((xg - 1j) * (hat * (E @ cp12))
                                - (xg + 1j) * (hat * (E @ cm12)))
    Fv = np.array([])
    def F_complex(z):
        s = mp.mpc(0.5 + z.imag, -z.real)
        val = 0.5 * s * (s - 1) * mp.pi ** (-s / 2) * mp.gamma(s / 2) \
            * mp.zeta(s)
        dval = mp.diff(lambda t: 0.5 * t * (t - 1) * mp.pi ** (-t / 2)
                       * mp.gamma(t / 2) * mp.zeta(t), s)
        return complex(z**2 * val / dval)
    Fv = np.array([F_complex(z) for z in xg])
    ratio = np.abs(W / Fv) / np.abs(R0_exact(A_NEW, np.pi, xg))
    b4 = {"median_M_over_R0_y0.3": float(np.median(ratio))}
    print(f"  B4 median |M/R0| at y0=0.3: {b4['median_M_over_R0_y0.3']:.3e}",
          flush=True)

    report.update({"B1_zeroflow": b1, "B2_gamma_locking": b2,
                   "B3_strip_winding": b3, "B4_envelope": b4})

    checks = {
        "VA1_backend_matches": va1 < 1e-6,
        "VA2_pipeline_a10": va2_q < 1e-6 and va2_z1 < 1e-6
        and va2_cnt == 0,
        "B3_no_new_complex_zeros": all(
            abs(b3[t][y]["excess"]) < 2.5 for t in b3 for y in b3[t]),
    }
    with CSV.open("w", encoding="utf-8", newline="") as f:
        f.write("series,a,j,value\n")
        for j, z in enumerate(zeros10):
            f.write(f"zero,10,{j+1},{z:.10f}\n")
        for j, z in enumerate(zeros12):
            f.write(f"zero,12,{j+1},{z:.10f}\n")
        for k, g in enumerate(gammas):
            f.write(f"gamma,-,{k+1},{g:.10f}\n")
    meta = {
        "experiment": "E-019",
        "task": "2.8(ii): backend v2 validation + zero-flow at a=12",
        "script": "research/numerics/scripts/E-019-zeroflow-a12.py",
        "parameters": {"a_val": A_VAL, "a_new": A_NEW, "h": H,
                       "lambda": LAM, "z_max": Z_MAX},
        "checks": checks,
        "report": report,
        "elapsed_seconds": round(time.time() - t0, 1),
        "timestamp_utc": datetime.now(timezone.utc).isoformat(),
    }
    JSON.write_text(json.dumps(meta, indent=2), encoding="utf-8")
    ok = all(checks.values())
    print(f"wrote {CSV.name}, {JSON.name} — "
          f"{'ALL CHECKS PASS' if ok else 'CHECK FAILURE: ' + str(checks)}")
    return 0 if ok else 1


if __name__ == "__main__":
    raise SystemExit(main())