E-014-cor16-diagnostics.py

{
  "experiment": "E-014",
  "task": "Q-S3: Corollary 1.6 diagnostics (Suzuki line)",
  "description": "theta-independence of zero count; v± localization; zero-flow map vs a; imaginary-axis least-squares read-off of theta(a), phi(a,z)",
  "script": "research/numerics/scripts/E-014-cor16-diagnostics.py",
  "parameters": {
    "a_sweep": [
      4.0,
      6.0,
      8.0,
      10.0
    ],
    "h": 0.01,
    "lambda": -1.0,
    "y_grid": [
      1.0,
      20.0,
      96
    ],
    "phi_deg": 4
  },
  "checks": {
    "V1_theta_count_interlacing": true,
    "V2_lattice_persists": true,
    "V3_xi_crosscheck_and_readoff": true,
    "V4_vmirror_symmetry": true
  },
  "key_values": {
    "P1_counts": [
      -76,
      -76,
      -76,
      -76,
      -76,
      -75,
      -75,
      -75,
      -75,
      -76,
      -76,
      -76,
      -76,
      -76,
      -76,
      -76
    ],
    "P1_spread": 1,
    "P2_localization": {
      "4.0": [
        0.8113352211750133,
        0.8042463400966671
      ],
      "6.0": [
        0.815502234769262,
        0.8096659632054103
      ],
      "8.0": [
        0.8230467697557222,
        0.8186214210851552
      ],
      "10.0": [
        0.8289899043479302,
        0.8251449076417641
      ]
    },
    "P3_first_zeros": {
      "4.0": [
        0.78678146,
        1.57356499,
        2.36035269,
        3.14714677,
        3.93394961
      ],
      "6.0": [
        0.52407766,
        1.04815542,
        1.57223347,
        2.09631224,
        2.62039153
      ],
      "8.0": [
        0.3929738,
        0.78594774,
        1.17892187,
        1.57189604,
        1.96487117
      ],
      "10.0": [
        0.31431482,
        0.62862966,
        0.94294449,
        1.2572596,
        1.57157429
      ]
    },
    "P4_readoff": {
      "4.0": {
        "best_theta": 3.141592653589793,
        "best_resid": 0.017714374721797224,
        "best_phase_agree": 1.0,
        "median_resid": 0.022342865357109994
      },
      "6.0": {
        "best_theta": 3.141592653589793,
        "best_resid": 0.021508567165581402,
        "best_phase_agree": 1.0,
        "median_resid": 0.0244687465471507
      },
      "8.0": {
        "best_theta": 3.141592653589793,
        "best_resid": 0.023225193635352376,
        "best_phase_agree": 1.0,
        "median_resid": 0.025389981806922614
      },
      "10.0": {
        "best_theta": 3.141592653589793,
        "best_resid": 0.024192319237001752,
        "best_phase_agree": 1.0,
        "median_resid": 0.025897228710069387
      }
    },
    "xi_flint_crosscheck_rel": 4.721803879534242e-17
  },
  "caveats": [
    "Galerkin discretization h=0.01 (E-010 h-convergence: z1 rel drift 4e-4 between h=0.01 and 0.005 at a=4)",
    "read-off on the imaginary axis tests only the modulus shape + phase of e^phi W vs F away from zeros; zero-density gap (F4 of D-020) means a good fit is NOT expected at large a — the residual trend is the observable"
  ],
  "elapsed_seconds": 20.4,
  "timestamp_utc": "2026-08-18T09:22:09.120726+00:00",
  "row_count": 81
}

"""E-014 — Q-S3 attack: Corollary 1.6 diagnostics on the Suzuki line.

Theory frame (D-023): Cor 1.6 demands e^{phi}W(a,theta;z) -> z^2 xi/xi'
uniformly on compacts. W has exponential type 2a with ALL REAL zeros
(Thm 1.5); at a <= 8 (E-010) the zeros form a maximal-density lattice
(spacing pi/a, i.e. ~2aT/pi zeros in [0,T]), while the target has
~T log T/2pi zeros. For the limit to hold, the zero density must COLLAPSE
at larger a — zeros stay real (Thm 1.5) so they must FLOW OUT to infinity
as a grows. Since theta only rotates the boundary phase (de Branges
interlacing), the count in [0,T] is theta-independent up to O(1); the
thinning, if it happens, lives in the a-dependence of the deficiency
vectors v± = T_a^{-1} e^{±x}. Four diagnostics:

  P1. theta-independence: zero count in [0, Z_TH] at a=6 for 16 thetas.
  P2. v± localization: edge-mass fraction of v± vs a (thinning channel:
      if v± concentrate at ±a, W's zeros can thin below maximal density).
  P3. zero-flow map: first 10 zeros of W at theta=0 for a in A_SWEEP;
      lattice law z_j ~ (pi/a)(1+rho)j and flow direction (inflow toward
      0 expected while the lattice persists; reversal = signal).
  P4. imaginary-axis read-off: F(iy) = -y^2 xi(1/2+y)/xi'(1/2+y) is real
      for real y. Fit log|W(a,theta;iy)| + p(y) ~ log|F(iy)| over
      y in [1, Y_MAX], p polynomial deg <= 4 (lstsq), theta on a 24-grid.
      Reports best theta(a) and residual vs a; phase-agreement score.

Checks:
  V1. P1: max count - min count <= 2 across theta (interlacing picture).
  V2. P3: lattice law persists: |z_1(a) - pi/a| / (pi/a) < 0.03 for all a.
  V3. P4: pipeline sanity — F(iy) cross-checked vs flint xi (rel 1e-8);
      read-off produces finite residuals with best <= median theta cost.
  V4. P2: v+ = mirror of v- (symmetry, rel 1e-6), machinery regression.

Outputs: research/numerics/data/E-014-cor16-diagnostics.csv (+ .json)
"""
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.linalg
import scipy.optimize
from numpy.polynomial.legendre import leggauss

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

T_MAXG = 20.0            # sieve to e^20 covers 2a <= 20 (a <= 10)
T_LERCH = 0.05
GL_N = 8
A_SWEEP = [4.0, 6.0, 8.0, 10.0]
H_MAIN = 0.01
LAM_SHIFT = -1.0
Z_TH = 40.0              # P1 window
N_GRID = 20000
N_ZEROS_TRACK = 10
THETAS_P1 = np.linspace(0, 2 * np.pi, 16, endpoint=False)
THETAS_P4 = np.linspace(0, 2 * np.pi, 24, endpoint=False)
Y_GRID = np.linspace(1.0, 20.0, 96)
PHI_DEG = 4


# ---------------------------------------------------------- g (as E-010)
def mangoldt_sieve_fast(n_max: int) -> np.ndarray:
    """Lambda(n)=log p for prime powers; vectorized Eratosthenes +
    vectorized prime-power assignment."""
    is_comp = np.zeros(n_max + 1, dtype=bool)
    for p in range(2, int(n_max ** 0.5) + 1):
        if not is_comp[p]:
            is_comp[p * p::p] = True
    lam_log = np.zeros(n_max + 1)
    primes = np.flatnonzero(~is_comp[2:]) + 2
    lp = np.log(primes)
    pk = primes.copy()
    while pk.size and pk[0] <= n_max:
        lam_log[pk] = lp
        pk = pk * primes[: pk.size]
        keep = pk <= n_max
        pk = pk[keep]
        lp = lp[keep]
    return lam_log


print("sieving von Mangoldt to e^20 ...", flush=True)
t_sv = time.time()
N_MAXG = int(np.exp(T_MAXG)) + 2
LAM_LOG = mangoldt_sieve_fast(N_MAXG)
NS = np.arange(2, N_MAXG + 1, dtype=np.float64)
LAM = LAM_LOG[2:] / np.sqrt(NS)
LOG_NS = np.log(NS)
C1 = np.cumsum(LAM)
C0 = np.cumsum(LAM * LOG_NS)
del LAM, LOG_NS, NS
print(f"  sieve done: {time.time()-t_sv:.0f}s", flush=True)

mp.mp.dps = 30
DIG_COEFF = float(-(mp.digamma(mp.mpf(1) / 4) - mp.log(mp.pi)) / 2)
PHI1 = float(mp.lerchphi(1, 2, mp.mpf(1) / 4))
KK = np.arange(0, 400)
DENOM = (KK + 0.25) ** 2


def prime_sum_vec(t):
    """sum_{n <= e^t} Lambda(n)/sqrt(n) * (t - log n), via C1/C0.
    #{2 <= n <= e^t} = ceil(e^t - eps) - 1 (searchsorted-right equivalent)."""
    m = np.floor(np.exp(t) + 1e-9).astype(np.int64) - 1
    out = np.zeros_like(t)
    pos = m > 0
    idx = m[pos] - 1
    out[pos] = C1[idx] * t[pos] - C0[idx]
    return out


def lerch_comb(t):
    out = np.empty_like(t)
    fast = t > T_LERCH
    if fast.any():
        z = np.exp(-2.0 * t[fast])
        phi = (np.power.outer(z, KK) / DENOM).sum(axis=1)
        out[fast] = -0.25 * (PHI1 - np.exp(-t[fast] / 2) * phi)
    if (~fast).any():
        for i in np.flatnonzero(~fast):
            ti = mp.mpf(float(t[i]))
            out[i] = float(-0.25 * (mp.lerchphi(1, 2, mp.mpf(1) / 4)
                                    - mp.exp(-ti / 2)
                                    * mp.lerchphi(mp.exp(-2 * ti), 2,
                                                  mp.mpf(1) / 4)))
    return out


def g_float(t):
    t = np.asarray(t, dtype=float)
    g = -4.0 * (np.exp(t / 2) + np.exp(-t / 2) - 2.0)
    g += prime_sum_vec(t)
    g += DIG_COEFF * t
    g += lerch_comb(t)
    g[t == 0.0] = 0.0
    return g


A_CONST = 0.5 * (np.log(2 * np.pi) + np.euler_gamma - 1.0)


def sigma_vec(t):
    t = np.asarray(t, dtype=float)
    with np.errstate(divide="ignore", invalid="ignore"):
        return np.where(t > 0, 0.5 * t * np.log(t) + A_CONST * t, 0.0)


def F_sing(t):
    if t <= 0:
        return 0.0
    return (A_CONST * t**3 / 6.0 + t**3 * np.log(t) / 12.0
            - 5.0 * t**3 / 72.0)


def S_sing(k, h):
    return (F_sing((k + 1) * h) - 2.0 * F_sing(k * h)
            + F_sing(abs(k - 1) * h))


XN, XW = leggauss(GL_N)


def S_full(k, h):
    """E-010's S with an adaptive cap: prime-point splitting only when the
    window holds <= 64 log-points (small t, where derivative jumps are
    relatively large); otherwise plain two-panel Gauss (g is C^0, jumps in
    derivative are O(Lambda/sqrt(n)) — negligible at large t)."""
    def quad_plain(lo, hi, peak):
        tot = 0.0
        for a, b in ((lo, peak), (peak, hi)):
            if b - a <= 0:
                continue
            mid, half = 0.5 * (a + b), 0.5 * (b - a)
            xs = mid + half * XN
            rho = g_float(xs) - sigma_vec(xs)
            tot += half * float(np.dot(XW, (h - np.abs(xs - peak)) * rho))
        return tot

    def 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:
            return quad_plain(lo, hi, peak)
        pts = 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})
        tot = 0.0
        for p0, p1 in zip(pts, pts[1:]):
            if p1 - p0 <= 0:
                continue
            mid, half = 0.5 * (p0 + p1), 0.5 * (p1 - p0)
            xs = mid + half * XN
            rho = g_float(xs) - sigma_vec(xs)
            tot += half * float(np.dot(XW, (h - np.abs(xs - peak)) * rho))
        return tot

    if k == 0:
        return S_sing(0, h) + 2.0 * quad(0.0, h, 0.0)
    return S_sing(k, h) + quad((k - 1) * h, (k + 1) * h, k * h)


def stiffness_first_col(n_elem, h):
    n_q = n_elem - 1
    s = np.array([S_full(k, h) for k in range(n_elem)])
    q = np.empty(n_q)
    q[0] = (2.0 * s[0] - 2.0 * s[1]) / h**2
    for k in range(1, n_q):
        q[k] = (2.0 * s[k] - s[k + 1] - s[k - 1]) / h**2
    return q


_BUILD = {}


def build(a, h):
    if (a, h) in _BUILD:
        return _BUILD[(a, h)]
    n_elem = int(round(2 * a / h))
    n = n_elem - 1
    q = stiffness_first_col(n_elem, h)
    Q = scipy.linalg.toeplitz(q)
    M = scipy.linalg.toeplitz(np.r_[2 * h / 3, h / 6,
                                    np.zeros(max(n - 2, 0))])
    ii = np.arange(1, n_elem)
    x = -a + ii * h
    xg, wg = leggauss(8)
    b = {}
    for sgn in (+1.0, -1.0):
        mid1 = x - h / 2
        xs1 = mid1[:, None] + (h / 2) * xg[None, :]
        phi1 = (xs1 - (x - h)[:, None]) / h
        mid2 = x + h / 2
        xs2 = mid2[:, None] + (h / 2) * xg[None, :]
        phi2 = ((x + h)[:, None] - xs2) / h
        b[sgn] = ((h / 2) * (np.exp(sgn * xs1) * phi1) @ wg
                  + (h / 2) * (np.exp(sgn * xs2) * phi2) @ wg)
    _BUILD[(a, h)] = (x, h, Q, M, b)
    return _BUILD[(a, h)]


def solve_v(Q, M, b_s, lam, n):
    T = Q - lam * M
    e = np.ones(n)
    B = np.zeros((n + 1, n + 1))
    B[:n, :n] = T
    B[:n, n] = e
    B[n, :n] = e
    sol = scipy.linalg.solve(B, np.r_[b_s, 0.0], assume_a="sym")
    return sol[:n], T


def hat_FT(z, h):
    return h * np.sinc(z * h / (2 * np.pi)) ** 2


def V_of_z(z, c, x, h):
    return hat_FT(z, h) * (np.exp(1j * np.outer(z, x)) @ c)


def W_complex(z, Vp, Vm, theta):
    return (z - 1j) * Vp + np.exp(1j * theta) * (z + 1j) * Vm


def zeros_theta0(a, h, lam, z_max, n_grid):
    """Zeros of W(a,0;z) in (0, z_max] — real gauge 2(z Re V + Im V)."""
    x, hh, Q, M, b = build(a, h)
    n = len(x)
    cp, T = solve_v(Q, M, b[+1.0], lam, n)
    zg = np.linspace(0.01, z_max, n_grid)
    Vp = V_of_z(zg, cp, x, hh)
    Wf = 2 * (zg * Vp.real + Vp.imag)

    def Wdirect(z):
        V = hat_FT(np.array([z]), hh)[0] * complex(np.exp(1j * z * x) @ cp)
        return 2 * (z * V.real + V.imag)

    idx = np.flatnonzero(np.sign(Wf[:-1]) * np.sign(Wf[1:]) < 0)
    roots = []
    for i in idx:
        try:
            roots.append(scipy.optimize.brentq(Wdirect, zg[i], zg[i + 1],
                                               xtol=1e-13, rtol=1e-13))
        except ValueError:
            pass
    return roots, cp, x, hh, Q, M, b, T


def count_zeros(a, h, lam, theta, z_max, n_grid):
    """Zero count in (0, z_max] via argument-principle-free |W| minima."""
    x, hh, Q, M, b = build(a, h)
    n = len(x)
    cp, T = solve_v(Q, M, b[+1.0], lam, n)
    cm, _ = solve_v(Q, M, b[-1.0], lam, n)
    np_ = float(cp @ T @ cp)
    nm = float(cm @ T @ cm)
    cm = cm * np.sqrt(np_ / nm)
    # argument-principle count on the rectangle [0.01, z_max] x [-eps, eps]
    # (all zeros real by Thm 1.5, so this counts real zeros in (0, z_max])
    eps = 0.02
    n_side = 4000
    xs = np.linspace(0.01, z_max, n_side)
    ys = np.linspace(-eps, eps, 41)
    contour = np.concatenate([
        xs + 1j * eps,          # top, left->right
        ys[::-1] * 1j + z_max,  # right, top->bottom
        xs[::-1] - 1j * eps,    # bottom, right->left
        ys * 1j + 0.01,         # left, bottom->top
    ])
    Vp = V_of_z(contour, cp, x, hh)
    Vm = V_of_z(contour, cm, x, hh)
    W = W_complex(contour, Vp, Vm, theta)
    darg = np.angle(W[1:] * np.conj(W[:-1])) + np.angle(W[0] * np.conj(W[-1]))
    count = int(round(float(darg.sum()) / (2 * np.pi)))
    return count, W


# ----------------------------------------------------------- F on i*y
def xi_real(s):
    """xi(s) for real s > 1 via mpmath (real-axis, safe)."""
    s = mp.mpf(s)
    return 0.5 * s * (s - 1) * mp.pi ** (-s / 2) * mp.gamma(s / 2) * mp.zeta(s)


mp.mp.dps = 40


def F_imag(y):
    """F(iy) = -y^2 xi(1/2+y)/xi'(1/2+y), real."""
    s = mp.mpf(0.5) + mp.mpf(y)
    xiv = xi_real(s)
    dx = mp.diff(xi_real, s)
    return float(-mp.mpf(y) ** 2 * xiv / dx)


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()
    rows = []

    # ---- V3 pre-check: flint cross-check of xi at two points
    import flint
    flint.ctx.prec = 200
    def xi_flint(s):
        a = flint.arb(str(s))
        return float(0.5 * a * (a - 1) * (a / 2).gamma()
                     * flint.arb.pi() ** (-a / 2) * a.zeta())
    xerr = float(max(abs(xi_real(2.5) - xi_flint(2.5)) / abs(xi_flint(2.5)),
                     abs(xi_real(10.5) - xi_flint(10.5))
                     / abs(xi_flint(10.5))))
    print(f"xi cross-check vs flint: rel err {xerr:.2e}", flush=True)

    Fv = np.array([F_imag(y) for y in Y_GRID])
    logF = np.log(np.abs(Fv))
    signF = np.sign(Fv)

    # ---- P1: theta-independence at a=6
    print("P1: theta sweep at a=6 ...", flush=True)
    counts = []
    for th in THETAS_P1:
        c, _ = count_zeros(6.0, H_MAIN, LAM_SHIFT, th, Z_TH, N_GRID)
        counts.append(c)
        print(f"  theta={th:.3f}  count={c}", flush=True)
    count_spread = max(counts) - min(counts)
    check_v1 = bool(count_spread <= 2)

    # ---- P2/P3/P4 over a-sweep
    flow = {}
    localization = {}
    readoff = {}
    sym_err = 0.0
    for a in A_SWEEP:
        print(f"a={a}: zeros + profiles + read-off ...", flush=True)
        roots, cp, x, hh, Q, M, b, T = zeros_theta0(a, H_MAIN, LAM_SHIFT,
                                                    max(Z_TH, 60.0), N_GRID)
        cm, _ = solve_v(Q, M, b[-1.0], LAM_SHIFT, len(x))
        np_ = float(cp @ T @ cp)
        nm = float(cm @ T @ cm)
        cm = cm * np.sqrt(np_ / nm)
        # V4: mirror symmetry v-(x) vs v+(-x)
        sym_err = max(sym_err,
                      float(np.max(np.abs(cm - cp[::-1]))
                            / max(1.0, float(np.max(np.abs(cp))))))
        # P2 localization: edge mass in T-norm and L2
        w_t = np.abs(T @ cp)
        edge = np.abs(x) > (a - 1.0)
        loc_t = float(np.sum(np.abs(cp[edge] * w_t[edge]))
                      / np.sum(np.abs(cp * w_t)))
        l2 = np.abs(cp) ** 2
        loc_l2 = float(l2[edge].sum() / l2.sum())
        localization[a] = (loc_t, loc_l2)
        # P3 flow
        flow[a] = roots[:N_ZEROS_TRACK]
        # P4 read-off
        ziy = 1j * Y_GRID
        Vp = V_of_z(ziy, cp, x, hh)
        Vm = V_of_z(ziy, cm, x, hh)
        best = None
        costs = []
        for th in THETAS_P4:
            W = W_complex(ziy, Vp, Vm, th)
            lw = np.log(np.abs(W))
            A = np.vander(Y_GRID, PHI_DEG + 1, increasing=True)
            coef, *_ = np.linalg.lstsq(A, logF - lw, rcond=None)
            resid = float(np.sqrt(np.mean((logF - lw - A @ coef) ** 2)))
            # phase agreement: arg(W) vs sign(F)
            phase = float(np.mean(np.cos(np.angle(W)) * signF > 0))
            costs.append((resid, th, phase))
        costs.sort()
        best = costs[0]
        med = float(np.median([c[0] for c in costs]))
        readoff[a] = {"best_theta": best[1], "best_resid": best[0],
                      "best_phase_agree": best[2], "median_resid": med}
        print(f"  loc(T)={loc_t:.4f} loc(L2)={loc_l2:.4f}  "
              f"z1={flow[a][0]:.4f} (pi/a={np.pi/a:.4f})  "
              f"read-off: theta*={best[1]:.3f} resid={best[0]:.4f} "
              f"(median {med:.4f}) phase={best[2]:.3f}", flush=True)

    # ---- checks
    lat_ok = all(abs(flow[a][0] - np.pi / a) / (np.pi / a) < 0.03
                 for a in A_SWEEP)
    check_v2 = bool(lat_ok)
    check_v3 = bool(xerr < 1e-8 and
                    all(readoff[a]["best_resid"]
                        <= readoff[a]["median_resid"] for a in A_SWEEP))
    check_v4 = bool(sym_err < 1e-6)
    checks = {
        "V1_theta_count_interlacing": check_v1,
        "V2_lattice_persists": check_v2,
        "V3_xi_crosscheck_and_readoff": check_v3,
        "V4_vmirror_symmetry": check_v4,
    }

    with CSV.open("w", encoding="utf-8", newline="") as f:
        f.write("part,a,theta_or_j,value\n")
        for th, c in zip(THETAS_P1, counts):
            f.write(f"P1_count,6.0,{th:.6f},{c}\n")
        for a in A_SWEEP:
            lt, l2 = localization[a]
            f.write(f"P2_loc_T,{a},-,{lt:.8f}\n")
            f.write(f"P2_loc_L2,{a},-,{l2:.8f}\n")
            for j, z in enumerate(flow[a], 1):
                f.write(f"P3_zero,{a},{j},{z:.10f}\n")
            ro = readoff[a]
            f.write(f"P4_best_theta,{a},-,{ro['best_theta']:.6f}\n")
            f.write(f"P4_best_resid,{a},-,{ro['best_resid']:.8f}\n")
            f.write(f"P4_median_resid,{a},-,{ro['median_resid']:.8f}\n")
            f.write(f"P4_phase_agree,{a},-,{ro['best_phase_agree']:.6f}\n")

    meta = {
        "experiment": "E-014",
        "task": "Q-S3: Corollary 1.6 diagnostics (Suzuki line)",
        "description": "theta-independence of zero count; v± localization; "
                       "zero-flow map vs a; imaginary-axis least-squares "
                       "read-off of theta(a), phi(a,z)",
        "script": "research/numerics/scripts/E-014-cor16-diagnostics.py",
        "parameters": {"a_sweep": A_SWEEP, "h": H_MAIN, "lambda": LAM_SHIFT,
                       "y_grid": [1.0, 20.0, 96], "phi_deg": PHI_DEG},
        "checks": checks,
        "key_values": {
            "P1_counts": counts, "P1_spread": count_spread,
            "P2_localization": {str(a): localization[a] for a in A_SWEEP},
            "P3_first_zeros": {str(a): [round(z, 8) for z in flow[a][:5]]
                               for a in A_SWEEP},
            "P4_readoff": {str(a): readoff[a] for a in A_SWEEP},
            "xi_flint_crosscheck_rel": float(xerr),
        },
        "caveats": [
            "Galerkin discretization h=0.01 (E-010 h-convergence: z1 rel "
            "drift 4e-4 between h=0.01 and 0.005 at a=4)",
            "read-off on the imaginary axis tests only the modulus shape "
            "+ phase of e^phi W vs F away from zeros; zero-density gap "
            "(F4 of D-020) means a good fit is NOT expected at large a — "
            "the residual trend is the observable",
        ],
        "elapsed_seconds": round(time.time() - t0, 1),
        "timestamp_utc": datetime.now(timezone.utc).isoformat(),
        "row_count": None,
    }
    meta["row_count"] = sum(1 for _ in CSV.open())
    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())