E-016-herglotz-fine-structure.py

{
  "experiment": "E-016",
  "task": "Q-S3c: Herglotz fine structure of W at theta=pi vs F",
  "description": "Poisson-smoothed zero density -(1/pi) Im W'/W(iy) at theta=pi (and theta=0 control) vs rho_F; analytic Galerkin z-derivative; zero-Poisson reconstruction cross-check",
  "script": "research/numerics/scripts/E-016-herglotz-fine-structure.py",
  "parameters": {
    "a_sweep": [
      4.0,
      6.0,
      8.0,
      10.0
    ],
    "h": 0.01,
    "lambda": -1.0,
    "y_grid": [
      0.5,
      20.0,
      391
    ],
    "z_max_zeros": 300.0
  },
  "checks": {
    "V1_plateau_density": true,
    "V2_rhoF_finitediff": true,
    "V3_poisson_reconstruction": false,
    "V4_vmirror_symmetry": true
  },
  "key_values": {
    "xi_flint_crosscheck_rel": 8.821833278620677e-17,
    "rhoF_finitediff_max_rel": 1.3322584092623682e-08,
    "metrics": {
      "4.0": {
        "n_zeros": 381,
        "z1": 1.1260734006983593,
        "plateau": 1.2713747663068067,
        "density_a_over_pi": 1.2732395447351628,
        "plateau_rel_dev": 0.0014645935527740487,
        "v3_rel": 0.32451368257591473,
        "kappa_shift": -0.21982980149584136,
        "d_pi": 1.2200734562385358,
        "d_zero": 1.2165204320831793
      },
      "6.0": {
        "n_zeros": 572,
        "z1": 0.7499016225074712,
        "plateau": 1.9083309455490012,
        "density_a_over_pi": 1.909859317102744,
        "plateau_rel_dev": 0.00080025347419903,
        "v3_rel": 0.2938167612479611,
        "kappa_shift": -0.39086861455767985,
        "d_pi": 1.855357135046592,
        "d_zero": 1.8529753375094085
      },
      "8.0": {
        "n_zeros": 763,
        "z1": 0.5621012628096893,
        "plateau": 2.5451204974114527,
        "density_a_over_pi": 2.5464790894703255,
        "plateau_rel_dev": 0.0005335178539225225,
        "v3_rel": 0.2792472337188209,
        "kappa_shift": -0.5618837715660804,
        "d_pi": 2.4914430957151077,
        "d_zero": 2.489482228965594
      },
      "10.0": {
        "n_zeros": 954,
        "z1": 0.44959700551214055,
        "plateau": 3.1817851862991344,
        "density_a_over_pi": 3.183098861837907,
        "plateau_rel_dev": 0.00041270334218086014,
        "v3_rel": 0.27086948306792724,
        "kappa_shift": -0.7328752211905654,
        "d_pi": 3.127620084707738,
        "d_zero": 3.126136732274815
      }
    },
    "vmirror_sym_err": 3.797860165317667e-14
  },
  "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)",
    "Poisson smoothing at height y averages over scale ~y: lattice oscillations at spacing pi/a are suppressed ~e^{-2ay}; the sensitive region is y <~ few (low zeros resolve individually)",
    "theta=pi vs theta=0 comparison is the sharp test of the H4 read-off (D-023): if theta=pi is correct, d_pi should lie systematically below d_zero as a grows"
  ],
  "elapsed_seconds": 32.7,
  "timestamp_utc": "2026-08-18T10:57:35.113504+00:00",
  "row_count": 3520
}

"""E-016 — Q-S3c: Herglotz (Poisson) fine structure of W at theta=pi vs F.

Theory (D-023 §4, item 3): for entire W with ALL REAL zeros z_j,
  -Im(W'/W)(iy) = -Im(kappa) + sum_j y/(y^2 + z_j^2)
is (pi times) the Poisson-smoothed zero density at height y. Individual
zeros within |z_j| <~ y leave distinguishable fingerprints, so this
diagnostic detects sub-lattice-scale zero redistribution BEFORE the
density collapses (Cor 1.6 requires the a->inf zeros to become the
gamma_j / xi'-zeros of F = z^2 xi(1/2-iz)/xi'(1/2-iz)).

For F on the imaginary axis (z = iy, s = 1/2+y real):
  rho_F(y) := -(1/pi) Im F'/F(iy)
           = (1/pi)[ 2/y + xi'/xi(1/2+y) - xi''/xi'(1/2+y) ].
Equivalently rho_F(y) = (1/pi) d/dy log|F(iy)|  (finite-diff check V2).

Diagnostics per a in {4,6,8,10} (h=0.01, lambda=-1, theta=pi; theta=0
control):
  D1. rho_W(y; a, pi) on y in [0.5, 20] via the ANALYTIC z-derivative of
      the Galerkin Fourier representation (no finite differences):
        V(z)  = hat(z) sum c_j e^{i z x_j}
        V'(z) = hat'(z) sum c_j e^{i z x_j} + hat(z) sum (i x_j) c_j e^{i z x_j}
        W(z)  = (z-i) Vp + e^{i theta} (z+i) Vm   (and z-derivative).
  D2. rho_F(y) analytic; cross-checked vs finite differences of
      log|F(iy)| (V2).
  D3. Zero-Poisson reconstruction: zeros of W(.,pi) located on the real
      axis (|W|-minima polished); rho from the Poisson sum must match D1
      up to the additive constant -Im(kappa)/pi (mean-matched over the
      window) — validates the derivative machinery (V3).
  D4. Metrics vs a: plateau value vs a/pi (density law); L2 distance
      d(a,theta) = ||rho_W - rho_F|| on [1,20] (theta=pi vs theta=0 —
      does theta=pi actually sharpen the approach?); oscillation
      envelope of rho_W at small y.

Checks:
  V1. plateau: |mean rho_W(y in [10,20]) - a/pi| / (a/pi) < 0.05.
  V2. rho_F analytic vs finite-diff of log|F|: max rel dev < 1e-6.
  V3. Poisson reconstruction from located zeros matches D1 (mean-matched
      L2 / rms < 0.02).
  V4. mirror symmetry v-(x) = v+(-x) rel err < 1e-6 (machinery regression).

Outputs: research/numerics/data/E-016-herglotz-fine-structure.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-016-herglotz-fine-structure.csv"
JSON = DATA / "E-016-herglotz-fine-structure.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
Y_GRID = np.linspace(0.5, 20.0, 391)
Y_PLATEAU = (10.0, 20.0)
Z_MAX_ZERO = 300.0
N_GRID = 40000


# ---------------------------------------------------------- g (as E-014)
def mangoldt_sieve_fast(n_max: int) -> np.ndarray:
    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):
    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):
    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


# ------------------------------------------- V(z), V'(z), W(z), W'(z)
def hat_FT(z, h):
    return h * np.sinc(z * h / (2 * np.pi)) ** 2


def hat_FT_d(z, h):
    """d/dz [h sinc(zh/2pi)^2], sinc(x)=sin(pi x)/(pi x)."""
    u = z * h / 2.0
    s = np.sinc(z * h / (2 * np.pi))          # sin(u)/u
    with np.errstate(divide="ignore", invalid="ignore"):
        ds = np.where(np.abs(u) > 1e-8,
                      (u * np.cos(u) - np.sin(u)) / u**2,
                      -u / 3.0)
    return h**2 * s * ds                      # d/dz[h s^2] = h^2 s ds/du


def V_and_dV(z, c, x, h):
    E = np.exp(1j * np.outer(z, x))
    S1 = E @ c
    S2 = (1j * x * E) @ c
    return hat_FT(z, h) * S1, hat_FT_d(z, h) * S1 + hat_FT(z, h) * S2


def W_and_dW(z, Vp, dVp, Vm, dVm, theta):
    eit = np.exp(1j * theta)
    W = (z - 1j) * Vp + eit * (z + 1j) * Vm
    dW = Vp + (z - 1j) * dVp + eit * (Vm + (z + 1j) * dVm)
    return W, dW


# --------------------------------------------------------------- F side
def xi_real(s):
    s = mp.mpf(s)
    if s == 1:                       # xi(1) = 1/2 (removable: zeta pole
        return mp.mpf(0.5)           # cancels the (s-1) factor)
    return 0.5 * s * (s - 1) * mp.pi ** (-s / 2) * mp.gamma(s / 2) * mp.zeta(s)


mp.mp.dps = 40


def rho_F_analytic(y):
    """rho_F(y) = (1/pi)[2/y + xi'/xi(1/2+y) - xi''/xi'(1/2+y)]."""
    s = mp.mpf(0.5) + mp.mpf(y)
    f = lambda t: xi_real(t)
    x0 = xi_real(s)
    x1 = mp.diff(f, s)
    x2 = mp.diff(f, s, 2)
    return float((2.0 / mp.mpf(y) + x1 / x0 - x2 / x1) / mp.pi)


def logabs_F_imag(y):
    s = mp.mpf(0.5) + mp.mpf(y)
    xiv = xi_real(s)
    dx = mp.diff(xi_real, s)
    return mp.log(abs(-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()

    # ---- F side
    print("F side: rho_F analytic + finite-diff cross-check ...", flush=True)
    # flint cross-check of xi_real (mpmath zeta is broken beyond ~16-19
    # digits in this env; E-014 verified 1e-8 at s=2.5,10.5 — extend to
    # the y-range endpoints used here)
    import flint
    flint.ctx.prec = 200
    def xi_flint(s):
        aa = flint.arb(str(s))
        return float(0.5 * aa * (aa - 1) * (aa / 2).gamma()
                     * flint.arb.pi() ** (-aa / 2) * aa.zeta())
    xi_err = float(max(abs(xi_real(s) - xi_flint(s)) / abs(xi_flint(s))
                       for s in (1.5, 10.5, 20.5)))  # s=1: arb zeta pole
                                                     # (xi(1)=1/2 exact)
    print(f"  xi vs flint max rel err: {xi_err:.2e}", flush=True)
    rhoF = np.array([rho_F_analytic(y) for y in Y_GRID])
    dlogF = np.array([
        float((logabs_F_imag(mp.mpf(y) + mp.mpf(1e-4))
               - logabs_F_imag(mp.mpf(y) - mp.mpf(1e-4))) / mp.mpf(2e-4))
        for y in Y_GRID]) / np.pi
    rel = np.abs(rhoF - dlogF) / np.maximum(np.abs(rhoF), 1e-12)
    v2_err = float(rel.max())
    print(f"  V2 finite-diff max rel dev: {v2_err:.2e}", flush=True)

    rows = []
    for y, r in zip(Y_GRID, rhoF):
        rows.append(f"F,{y:.6f},{r:.10e}")

    sym_err = 0.0
    metrics = {}
    for a in A_SWEEP:
        print(f"a={a}: solve v±, Herglotz densities ...", flush=True)
        x, hh, Q, M, b = build(a, H_MAIN)
        n = len(x)
        cp, T = solve_v(Q, M, b[+1.0], LAM_SHIFT, n)
        cm, _ = solve_v(Q, M, b[-1.0], LAM_SHIFT, n)
        np_ = float(cp @ T @ cp)
        nm = float(cm @ T @ cm)
        cm = cm * np.sqrt(np_ / nm)
        sym_err = max(sym_err,
                      float(np.max(np.abs(cm - cp[::-1]))
                            / max(1.0, float(np.max(np.abs(cp))))))

        ziy = 1j * Y_GRID
        Vp, dVp = V_and_dV(ziy, cp, x, hh)
        Vm, dVm = V_and_dV(ziy, cm, x, hh)

        out = {}
        for th, tag in ((np.pi, "pi"), (0.0, "zero")):
            W, dW = W_and_dW(ziy, Vp, dVp, Vm, dVm, th)
            rhoW = -np.imag(dW / W) / np.pi
            out[tag] = rhoW
            for y, r in zip(Y_GRID, rhoW):
                rows.append(f"W_th{tag},{a},{y:.6f},{r:.10e}")

        # ---- zeros of W(., pi) on the real axis (for V3)
        # c = sqrt(np_/nm) = 1 exactly (mirror symmetry), so for real z,
        # Vm = conj(Vp) and W(z) = 2i*g(z) with the REAL gauge
        # g(z) = z Im Vp - Re Vp. Zeros = sign changes of g -> brentq.
        assert abs(np_ / nm - 1.0) < 1e-12, "norm asymmetry breaks gauge"
        zg = np.linspace(1e-3, Z_MAX_ZERO, N_GRID)
        Vpr, _ = V_and_dV(zg, cp, x, hh)
        g = zg * Vpr.imag - Vpr.real

        def g_direct(zz):
            V = hat_FT(np.array([zz]), hh)[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
        zeros = np.array(zeros)

        # V3: Poisson reconstruction (+- symmetric located zeros to
        # Z_MAX_ZERO + Euler-Maclaurin lattice tail beyond, spacing pi/a,
        # phase matched to last located zero; mean-matched vs analytic).
        d_lat = np.pi / a
        z_last = zeros[-1] + d_lat / 2.0
        recon = np.empty_like(Y_GRID)
        for k, y in enumerate(Y_GRID):
            tail = (y / d_lat) * (np.pi / 2.0 - np.arctan(z_last / y))
            recon[k] = (2.0 * (np.sum(y / (y**2 + zeros**2)) + tail)) / np.pi
        rhoW = out["pi"]
        shift = float(np.mean(rhoW - recon))
        rms = float(np.sqrt(np.mean(rhoW**2)))
        v3_rel = float(np.sqrt(np.mean((rhoW - recon - shift)**2)) / rms)

        # D4 metrics
        pm = (Y_GRID >= Y_PLATEAU[0]) & (Y_GRID <= Y_PLATEAU[1])
        plateau = float(np.mean(rhoW[pm]))
        dens = a / np.pi
        sel = Y_GRID >= 1.0
        d_pi = float(np.sqrt(np.mean((out["pi"][sel] - rhoF[sel])**2)))
        d_zero = float(np.sqrt(np.mean((out["zero"][sel] - rhoF[sel])**2)))
        metrics[a] = {
            "n_zeros": len(zeros), "z1": float(zeros[0]) if len(zeros) else None,
            "plateau": plateau, "density_a_over_pi": dens,
            "plateau_rel_dev": abs(plateau - dens) / dens,
            "v3_rel": v3_rel, "kappa_shift": shift,
            "d_pi": d_pi, "d_zero": d_zero,
        }
        print(f"  zeros={len(zeros)} z1={metrics[a]['z1']:.4f} "
              f"plateau={plateau:.4f} (a/pi={dens:.4f}) "
              f"V3 rel={v3_rel:.2e} d_pi={d_pi:.4f} d_zero={d_zero:.4f}",
              flush=True)

    checks = {
        "V1_plateau_density": bool(all(
            metrics[a]["plateau_rel_dev"] < 0.05 for a in A_SWEEP)),
        "V2_rhoF_finitediff": bool(v2_err < 1e-6),
        "V3_poisson_reconstruction": bool(all(
            metrics[a]["v3_rel"] < 0.02 for a in A_SWEEP)),
        "V4_vmirror_symmetry": bool(sym_err < 1e-6),
    }

    with CSV.open("w", encoding="utf-8", newline="") as f:
        f.write("series,a,y,rho\n")
        for r in rows:
            parts = r.split(",")
            if parts[0] == "F":
                f.write(f"F,-,{parts[1]},{parts[2]}\n")
            else:
                f.write(f"{parts[0]},{parts[1]},{parts[2]},{parts[3]}\n")

    meta = {
        "experiment": "E-016",
        "task": "Q-S3c: Herglotz fine structure of W at theta=pi vs F",
        "description": "Poisson-smoothed zero density -(1/pi) Im W'/W(iy) "
                       "at theta=pi (and theta=0 control) vs rho_F; "
                       "analytic Galerkin z-derivative; zero-Poisson "
                       "reconstruction cross-check",
        "script": "research/numerics/scripts/E-016-herglotz-fine-structure.py",
        "parameters": {"a_sweep": A_SWEEP, "h": H_MAIN, "lambda": LAM_SHIFT,
                       "y_grid": [0.5, 20.0, 391], "z_max_zeros": Z_MAX_ZERO},
        "checks": checks,
        "key_values": {
            "xi_flint_crosscheck_rel": xi_err,
            "rhoF_finitediff_max_rel": v2_err,
            "metrics": {str(a): metrics[a] for a in A_SWEEP},
            "vmirror_sym_err": sym_err,
        },
        "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)",
            "Poisson smoothing at height y averages over scale ~y: lattice "
            "oscillations at spacing pi/a are suppressed ~e^{-2ay}; the "
            "sensitive region is y <~ few (low zeros resolve individually)",
            "theta=pi vs theta=0 comparison is the sharp test of the H4 "
            "read-off (D-023): if theta=pi is correct, d_pi should lie "
            "systematically below d_zero as a grows",
        ],
        "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())