E-021-sinf-sf-curve.py

{
  "experiment": "E-021",
  "task": "3.7: exact S_inf(n) + |S_f(n)| curve + D-036 constants",
  "script": "research/numerics/scripts/E-021-sinf-sf-curve.py",
  "parameters": {
    "n_max": 5000,
    "t_zeros": 5000.0
  },
  "checks": {
    "V1_lambda_vs_E015": true,
    "V2_R_bound_holds": true,
    "V3_Sinf_residual_bounded": true
  },
  "key_values": {
    "A_use": 0.9,
    "B0_use": 1.0,
    "R_max_sample": 4.7507714994958565,
    "C0_samples": {
      "1e2": 6.234438640727731,
      "1e4": 8.840233256550178,
      "1e6": 11.194201354663145
    },
    "eta_0": "-0.5772156649015328606065120900824024310",
    "Sinf_residual_fits": [
      [
        1000,
        -0.2500416666625
      ],
      [
        2000,
        -0.2500208333328125
      ],
      [
        3000,
        -0.25001388888873455
      ],
      [
        4000,
        -0.25001041666660156
      ],
      [
        5000,
        -0.2500083333333
      ]
    ],
    "V1_worst_rel": 4.79820218420929e-25,
    "S_f_max_abs": "13.5229134371",
    "S_f_max_abs_n": 3298
  },
  "elapsed_seconds": 2139.3,
  "timestamp_utc": "2026-08-19T06:48:55.463113+00:00"
}

"""E-021 — task 3.7: exact S_inf(n) to n=1e4, |S_f(n)| curve, tightened
D-036 constants. Feeds D-037 (finite-height Li theorem).

Components:
  P1. D-036 constants: C0(t) from Trudgian; |R(s)| = |zeta'/zeta(s) -
      sum_{|gamma-t|<1} 1/(s-rho)| sampled at ~30 contour points (flint
      acb zeta + fd derivative, zeros via mp.zetazero) -> fit A, B0.
  P2. eta_j — ALL FLINT (session 006 cont.10 root-cause: mpmath zeta,
      Hurwitz included, silently loses accuracy beyond ~58 digits on
      this stack [mpmath 1.3 + Py 3.14] -> eta_j errors ~1e-58 ->
      S_f binomial cancellation explodes ~(4/3)^n; V1 failed 1.4e143.
      The zero-side sigma route also has eta_0 wrong by 2.4e-4 from
      sigma_1 truncation at T=5000 — both defects removed here).
      j <= 23: Stieltjes series division (E-011 machinery, validated):
        -zeta'/zeta(1+u) = (1/u) N(u)/D(u), eta_k = [u^{k+1}] N/D.
      j >= 24: Lagarias (4.10) with flint arb zeta:
        eta_j = (-1)^{j+1} [(1-2^{-(j+1)}) zeta(j+1) - 1 + sigma_{j+1}]
        with sigma_m = (-1)^{m-1}(tau_{m-1} - eta_{m-1}) + 1 EXACT
        (arithmetic side) for m <= 24, and sigma_m = 2 Re sum
        rho^{-m} over the first 100 zeros (mp.zetazero, dps 100)
        for 25 <= m <= M_SIG; omitted beyond (boundary error
        C(n,M_SIG)*14^{-M_SIG} < 1e-40 by construction).
      NOTES (hard-won, session 006 cont.10):
      * sigma_m may NOT be omitted at any m <= M_SIG — the truncation
        boundary error ~ C(n,M)*14^{-M} PEAKS at M = n/15 with size
        ~(15/14)^n: with M=200, S_f exploded as 1e133 at n=5000
        (measured). M_SIG ~ 0.22*n_max.
      * the m=24 omission alone moved S_f(200) by ~0.4 (measured).
      * float64 zeros for sigma_m are NOT enough either (S_f(300) off
        by 7.8e-8 — measured); need mp.zetazero precision. Zeros
        k > 100 contribute < (14.13/236.5)^25 ~ 1e-31 relative for
        m >= 25, so the first 100 at dps 100 are effectively exact.
      * GAMMA_E must be a string (bare float -> 53-bit, 4.94e-21).
  P3. S_inf(n) EXACT for n <= N_MAX (Lagarias (4.11)):
        S_inf(n) = sum_{j=1}^n (-1)^j C(n,j) (1-2^{-j}) zeta*(j),
        zeta*(1) = log(4pi)+gamma, zeta*(j) = zeta(j) (flint.arb).
      Cancellation ~2^n -> dps = ceil(0.31*n) + 60.
  P4. S_f(n) = sum_{j=1}^n C(n,j) eta_{j-1} (plain binomial — the
      E-011/Lagarias-(4.13) convention after the (-1) absorbed into
      the series-eta definition; VALIDATED: S_f(200) = -2.14037668
      matches S_inf(200)+1-lambda_200(E-015)); lambda_n =
      S_inf - S_f + 1.
  P5. VALIDATION: lambda_n (n <= 3000) vs E-015 certified CSV (rel dev);
      eta_0 vs E-011 value; S_inf fit: S_inf - [(1/2)n log n + C1 n]
      -> c0 (+ c1/n...), C1 = (gamma-1-log(2pi))/2.

Checks:
  V1. lambda_n matches E-015 for n <= 3000 (rel dev < 1e-50 at dps used).
  V2. |R(s)| <= A log(t+2) + B0 at all sample points with the FITTED
      A, B0 (and fitted values reported vs D-036's 2.0/2.0).
  V3. S_inf fit residual -> constant c0 (|residual - c0| -> 0 trend).

Outputs: research/numerics/data/E-021-sinf-sf-curve.csv (+ .json)
"""
import json
import sys
import time
from datetime import datetime, timezone
from pathlib import Path

import flint
import mpmath as mp
import numpy as np

REPO = Path(__file__).resolve().parents[3]
DATA = REPO / "research" / "numerics" / "data"
CSV = DATA / "E-021-sinf-sf-curve.csv"
JSON = DATA / "E-021-sinf-sf-curve.json"
E015 = DATA / "E-015-li-coefficients-n3000.csv"

N_MAX = 5000            # 10000 costs ~2h (row cost ~ n^2 * dps);
                        # cached zeta/zeros make a later 1e4 rerun cheap
T_ZEROS = 5000.0        # zeros up to this height for sigma sums
N_ZMAX = 8000           # mp.zetazero count cap (N(5e3) ~ 7.9e3)
DPS_BIG = 3300
# NOTE: must be a STRING — a bare float literal collapses to 53-bit
# (measured error 4.94e-21 -> sigma_1/zeta*(1) wrong at digit 21).
GAMMA_E = "0.57721566490153286060651209008240243104215933593992337"


# ---------------------------------------------------------------- N(T)
def N_count(T):
    F = (T / (2 * np.pi)) * np.log(T / (2 * np.pi)) - T / (2 * np.pi) + 7 / 8
    return F


def S_bound(T):
    return 0.137 * np.log(T) + 0.443 * np.log(np.log(T)) + 1.588


def C0(t):
    """Upper bound for zeros in [t, t+1] (t >= 10)."""
    return (N_count(t + 1) - N_count(t)) + S_bound(t + 1) + S_bound(t)


# ------------------------------------------------------------- zeros
def compute_zeros():
    cache = DATA / "E-021-zeros-cache.npy"
    n_est = int(N_count(T_ZEROS) + S_bound(T_ZEROS)) + 100
    n_est = min(n_est, N_ZMAX)
    if cache.exists():
        zs = np.load(cache)
        if len(zs) >= n_est:
            print(f"  loaded {n_est} zeros from cache", flush=True)
            return zs[:n_est]
    print(f"  computing {n_est} zeros via mp.zetazero ...", flush=True)
    mp.mp.dps = 30
    zs = []
    for k in range(1, n_est + 1):
        zs.append(mp.zetazero(k).imag)
        if k % 1000 == 0:
            print(f"    {k}/{n_est}", flush=True)
    out = np.array([float(z) for z in zs])
    np.save(cache, out)
    return out


# ------------------------------------------------------- zeta strings
def get_zeta_strings(n_max):
    """zeta(j) decimal strings at full precision (flint arb), cached.
    full_digits = 0.31*n_max + 100 covers the ~2^n cancellation in
    BOTH S_inf and the (1-2^{-s})zeta(s)-1 eta formula (absolute
    error ~1e-1650 uniformly)."""
    full_digits = int(0.31 * n_max) + 100
    zcache = DATA / "E-021-zeta-cache.json"
    if zcache.exists():
        zeta_str = {int(k): v for k, v in
                    json.loads(zcache.read_text()).items()}
        print("  zeta strings loaded from cache", flush=True)
        # patch legacy zeta*(1) (built with double-precision GAMMA_E)
        mp.mp.dps = full_digits + 50
        z1 = mp.nstr(mp.log(4 * mp.pi) + mp.mpf(GAMMA_E), full_digits)
        if zeta_str.get(1) != z1:
            zeta_str[1] = z1
            zcache.write_text(json.dumps(zeta_str))
            print("  patched zeta*(1) in cache (GAMMA_E float fix)",
                  flush=True)
        return zeta_str, full_digits
    flint.ctx.prec = int(full_digits * 3.33) + 200
    mp.mp.dps = full_digits + 50
    zeta_str = {1: mp.nstr(mp.log(4 * mp.pi) + mp.mpf(GAMMA_E),
                           full_digits)}
    for j in range(2, n_max + 1):
        zeta_str[j] = flint.arb(j).zeta().str(full_digits,
                                              radius=False)
        if j % 2000 == 0:
            print(f"    zeta({j}) done", flush=True)
    zcache.write_text(json.dumps(zeta_str))
    return zeta_str, full_digits


# ---------------------------------------------------------------- eta
# K_SPLIT: eta_j for j < K_SPLIT from Stieltjes (exact, arithmetic);
# sigma_m exact for m <= K_SPLIT. Must be ~80: with K_SPLIT=24 the
# zero-tail (k>100) of sigma_25 entered at term j=24 with C(n,24)
# amplification — S_f(3000) was off by 1.06e3 growing ~n^24 (measured).
# With K_SPLIT=80 the zero-tail exposure (m>=81, k>100) is < 3e-11.
K_SPLIT = 80


def compute_eta(zeta_str, full_digits):
    """eta_j for j = 0..N_MAX-1, ALL flint-based (see header P2).
    Returns dict j -> decimal string at full_digits."""
    print(f"  eta_0..{K_SPLIT - 1}: flint Stieltjes + series division "
          "(E-011) ...", flush=True)
    flint.ctx.prec = int((full_digits + 100) * 3.33) + 200
    mp.mp.dps = full_digits + 100   # BEFORE mp.mpf(string) parses below!
    gam_s = [mp.mpf(flint.acb.stieltjes(k).real.str(
        full_digits + 100, radius=False)) for k in range(K_SPLIT + 2)]
    N_ = [mp.mpf(0)] * (K_SPLIT + 2)
    D_ = [mp.mpf(0)] * (K_SPLIT + 2)
    N_[0] = D_[0] = 1
    for m in range(0, K_SPLIT + 1):
        D_[m + 1] = (-1) ** m * gam_s[m] / mp.factorial(m)
    for m in range(1, K_SPLIT + 1):
        N_[m + 1] = -(-1) ** m * gam_s[m] / mp.factorial(m - 1)
    Ev = [mp.mpf(0)] * (K_SPLIT + 2)
    for k in range(0, K_SPLIT + 2):
        s = N_[k]
        for i in range(1, k + 1):
            s -= D_[i] * Ev[k - i]
        Ev[k] = s
    eta_str = {j: mp.nstr(Ev[j + 1], full_digits)
               for j in range(K_SPLIT)}
    print(f"    eta_0 = {mp.nstr(Ev[1], 15)} (should be -gamma)",
          flush=True)
    print(f"    eta_1 = {mp.nstr(Ev[2], 15)}", flush=True)

    # ---- sigma_m: exact (arithmetic) for m <= 24; zeros for 25..200
    print("  sigma_m: m<=24 exact (tau-eta), 25..200 from zeros ...",
          flush=True)
    mp.mp.dps = full_digits + 50
    tau = {0: mp.mpf(-0.5) * (mp.log(4 * mp.pi) + mp.mpf(GAMMA_E))}
    for j in range(1, K_SPLIT):
        z = mp.mpf(zeta_str[j + 1])
        tau[j] = (-1) ** (j + 1) * (1 - mp.mpf(2) ** (-(j + 1))) * z
    sigma = {}
    for m in range(1, K_SPLIT + 1):
        sigma[m] = ((-1) ** (m - 1)
                    * (tau[m - 1] - mp.mpf(eta_str[m - 1])) + 1)
    mp.mp.dps = 60
    sig1_known = 1 + mp.euler / 2 - mp.log(4 * mp.pi) / 2
    print(f"    sigma_1 = {mp.nstr(sigma[1], 20)} "
          f"(closed form {mp.nstr(sig1_known, 20)})", flush=True)
    # sigma_m, m = 25..M_SIG: first 100 zeros via mp.zetazero, dps 100.
    # M_SIG must be adaptive: truncating sigma at M leaves a boundary
    # error ~ C(n,M)*14.13^{-M} which PEAKS at M = n/15 (~(15/14)^n);
    # with M=200 the run exploded as 1e133 at n=5000 (measured). Need
    # C(n,M)*14^{-M} < 1e-40 on the decreasing side (M > n/15).
    # float64 zeros are NOT enough (S_f(300) off 7.8e-8 — measured);
    # zeros k>100 are suppressed by (14.13/236.5)^25 ~ 1e-31 relative.
    n_zeta = max(zeta_str)
    M_SIG = n_zeta
    for M in range(K_SPLIT + 1, n_zeta + 1):
        if M > n_zeta / 15:
            logC = (mp.loggamma(n_zeta + 1) - mp.loggamma(M + 1)
                    - mp.loggamma(n_zeta - M + 1))
            if float(logC - (M - 1) * mp.log(14.13)
                     + mp.log(0.05)) / np.log(10) < -40:
                M_SIG = M
                break
    print(f"    M_SIG = {M_SIG} (C(n,M)*14^{{-M}} < 1e-40)", flush=True)
    # sigma_m precision matters: dps-100 sums left eps_m ~ 14^{-m}e-100,
    # whose C(n,j)-amplified impact (peak term C(n,n/15)*14^{-n/15} ~
    # 1e105 at n=3000) showed up as S_f error ~1e3 (measured). Zeros
    # at 400 digits + sums at full_digits+50 kill it.
    # dps BEFORE any mp.mpf(string) parse! (Third instance of this bug
    # class: g100 parsed at leftover dps 60 corrupted sigma_m at
    # ~m*1e-60*14^{-m} -> S_f(3000) off by 1.7e7, growth ~n^80.)
    mp.mp.dps = full_digits + 50
    z100cache = DATA / "E-021-zeros100-hi-cache.json"
    if z100cache.exists():
        g100 = [mp.mpf(v) for v in
                json.loads(z100cache.read_text())]
    else:
        mp.mp.dps = 450
        g100 = [mp.zetazero(k).imag for k in range(1, 101)]
        z100cache.write_text(json.dumps([mp.nstr(g, 400)
                                         for g in g100]))
        mp.mp.dps = full_digits + 50
    rho = [mp.mpc(mp.mpf("0.5"), g) for g in g100]
    for m in range(K_SPLIT + 1, M_SIG + 1):
        sigma[m] = 2 * mp.re(sum(r ** (-m) for r in rho))
    mp.mp.dps = 60
    print(f"    sigma_2 = {mp.nstr(sigma[2], 20)}", flush=True)
    print(f"    sigma_25 (zeros) = {mp.nstr(sigma[25], 6)}; "
          f"sigma_{M_SIG} = {mp.nstr(sigma[M_SIG], 6)}", flush=True)

    print("  eta_j (j>=24): (-1)^{j+1}[(1-2^{-(j+1)})zeta(j+1) - 1 "
          "+ sigma_{j+1}] ...", flush=True)
    mp.mp.dps = full_digits + 50
    for j in range(K_SPLIT, n_zeta):
        z = mp.mpf(zeta_str[j + 1])
        sig = sigma[j + 1] if j + 1 <= M_SIG else mp.mpf(0)
        eta = (-1) ** (j + 1) * ((1 - mp.mpf(2) ** (-(j + 1))) * z - 1
                                 + sig)
        eta_str[j] = mp.nstr(eta, full_digits)
    # stitch check: eta_{K_SPLIT-1} via Stieltjes vs via full formula
    mp.mp.dps = full_digits + 50
    eta23_formula = ((-1) ** K_SPLIT
                     * ((1 - mp.mpf(2) ** (-K_SPLIT))
                        * mp.mpf(zeta_str[K_SPLIT])
                        - 1 + sigma[K_SPLIT]))
    stitch = abs(eta23_formula - mp.mpf(eta_str[K_SPLIT - 1]))
    mp.mp.dps = 60
    print(f"    stitch |eta_{K_SPLIT - 1} stieltjes - formula| = "
          f"{mp.nstr(stitch, 3)}", flush=True)
    return eta_str


# ----------------------------------------------------------- S_inf, Sf
def compute_S_inf(n_max, zeta_str, full_digits):
    """Exact S_inf(n) via zeta*(j) (flint arb), dps scaled by n.
    Cancellation ~2^n: every zeta(j) stored at FULL string precision
    (the Sigma(-1)^j C(n,j)*1 = 0 piece must cancel exactly)."""
    print("  parsing zeta strings once at full precision ...", flush=True)
    mp.mp.dps = full_digits + 50
    zeta_mpf = {j: mp.mpf(v) for j, v in zeta_str.items()}
    results = {}
    for n in range(1, n_max + 1):
        mp.mp.dps = int(0.31 * n) + 80
        s = mp.mpf(0)
        c = mp.mpf(1)
        for j in range(1, n + 1):
            c = c * (n - j + 1) / j          # C(n,j) by recurrence
            s += ((-1) ** j) * c * (1 - mp.mpf(2) ** (-j)) * zeta_mpf[j]
        results[n] = +s
        if n % 500 == 0:
            print(f"    S_inf row {n}/{n_max}", flush=True)
    return results


def compute_Sf(n_max, eta_str):
    """S_f(n) = sum_j C(n,j) eta_{j-1} — plain binomial (E-011/Lagarias
    (4.13) convention with series-eta; validated at n=200 vs E-015:
    S_f(200) = -2.14037668 = S_inf(200)+1-lambda_200). Requires eta_j
    absolute accuracy ~1e-1650 at n=5000 (terms C(n,j)*eta up to
    ~1e312 cancel to O(1)) — flint-only computation, see P2."""
    full_digits = int(0.31 * n_max) + 100
    mp.mp.dps = full_digits + 50
    eta_mpf = {j: mp.mpf(v) for j, v in eta_str.items()}
    results = {}
    for n in range(1, n_max + 1):
        mp.mp.dps = int(0.31 * n) + 80
        s = mp.mpf(0)
        c = mp.mpf(1)
        for j in range(1, n + 1):
            c = c * (n - j + 1) / j
            s += c * eta_mpf[j - 1]
        results[n] = +s
        if n % 500 == 0:
            print(f"    S_f row {n}/{n_max}", flush=True)
    return results


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

    # ---- P1: constants
    print("P1: D-036 constants ...", flush=True)
    print(f"  C0(100) = {C0(100):.3f}, C0(1e4) = {C0(1e4):.3f}, "
          f"C0(1e6) = {C0(1e6):.3f}", flush=True)
    gam_sample = compute_zeros()

    def xi_mp(s):
        return 0.5 * s * (s - 1) * mp.pi ** (-s / 2) * mp.gamma(s / 2) \
            * mp.zeta(s)

    def zeta_logder(s):
        """zeta'/zeta via xi'/xi minus explicit terms (mpmath dps 30;
        complex-zeta OK to ~1e-10 per E-018 C2)."""
        sm = mp.mpc(s.real, s.imag)
        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))

    R_samples = []
    mp.mp.dps = 30
    for sig in (-2.5, 0.5, 3.0):
        for t in (10.0, 50.0, 100.0, 500.0, 1000.0):
            s = complex(sig, t)
            near = [g for g in gam_sample if abs(g - t) < 1.0]
            psum = sum(1.0 / (s - complex(0.5, g)) for g in near)
            R_samples.append((sig, t, abs(zeta_logder(s) - psum)))
    # fit |R| <= A log(t+2) + B0
    A_fit = max(v / np.log(t + 2) for _, t, v in R_samples)
    A_use = max(0.5, float(np.ceil(A_fit * 1.2 * 10) / 10))
    B0_need = max(v - A_use * np.log(t + 2) for _, t, v in R_samples)
    B0_use = max(1.0, float(np.ceil(B0_need * 1.3 * 10) / 10))
    print(f"  |R| max sample = {max(v for *_ , v in R_samples):.2f}; "
          f"A_use={A_use}, B0_use={B0_use}", flush=True)

    # ---- P2: zeta strings + eta (all flint)
    print("P2: zeta strings + eta_j (flint) ...", flush=True)
    zeta_str, full_digits = get_zeta_strings(N_MAX)
    eta_str = compute_eta(zeta_str, full_digits)

    # ---- P3/P4: S_inf, S_f
    print("P3: S_inf exact ...", flush=True)
    S_inf = compute_S_inf(N_MAX, zeta_str, full_digits)
    print("P4: S_f ...", flush=True)
    S_f = compute_Sf(N_MAX, eta_str)

    # ---- P5: validation + fit
    print("P5: validation vs E-015 ...", flush=True)
    e015 = {}
    with E015.open() as f:
        header = f.readline()
        for line in f:
            parts = line.strip().split(",")
            if len(parts) >= 2 and parts[0].isdigit():
                e015[int(parts[0])] = mp.mpf(parts[1])
    worst = 0.0
    mp.mp.dps = 100
    for n in sorted(e015):
        lam = S_inf[n] - S_f[n] + 1
        dev = abs(lam - e015[n]) / abs(e015[n])
        worst = max(worst, float(dev))
    print(f"  V1 worst rel dev vs E-015 (n<=3000): {worst:.2e}", flush=True)

    C1 = (mp.mpf(GAMMA_E) - 1 - mp.log(2 * mp.pi)) / 2
    fits = []
    for n in (1000, 2000, 3000, 4000, 5000, 6000, 8000, 10000):
        if n > N_MAX:
            continue
        mp.mp.dps = 100
        resid = S_inf[n] - (mp.mpf(n) / 2) * mp.log(n) - C1 * n
        fits.append((n, float(resid)))
    print(f"  S_inf residual (should -> c0): {fits}", flush=True)

    # ---- write CSV
    with CSV.open("w", encoding="utf-8", newline="") as f:
        f.write("n,S_inf,S_f,lambda_n,abs_Sf\n")
        for n in range(1, N_MAX + 1):
            lam = S_inf[n] - S_f[n] + 1
            f.write(f"{n},{mp.nstr(S_inf[n],12)},{mp.nstr(S_f[n],12)},"
                    f"{mp.nstr(lam,12)},{mp.nstr(abs(S_f[n]),12)}\n")

    checks = {
        # V1 floor is ~1e-24 (E-015 CSV holds 26 digits); 1e-15 still
        # kills any macroscopic (convention/precision) bug — the bugs
        # found in cont.10 produced 1e-8..1e143 deviations.
        "V1_lambda_vs_E015": bool(worst < 1e-15),
        "V2_R_bound_holds": bool(all(
            v <= A_use * np.log(t + 2) + B0_use for _, t, v in R_samples)),
        "V3_Sinf_residual_bounded": bool(
            abs(fits[-1][1] - fits[-2][1]) < abs(fits[0][1] - fits[1][1])),
    }
    meta = {
        "experiment": "E-021",
        "task": "3.7: exact S_inf(n) + |S_f(n)| curve + D-036 constants",
        "script": "research/numerics/scripts/E-021-sinf-sf-curve.py",
        "parameters": {"n_max": N_MAX, "t_zeros": T_ZEROS},
        "checks": checks,
        "key_values": {
            "A_use": A_use, "B0_use": B0_use,
            "R_max_sample": max(v for *_, v in R_samples),
            "C0_samples": {"1e2": C0(100), "1e4": C0(1e4), "1e6": C0(1e6)},
            "eta_0": eta_str[0][:40],
            "Sinf_residual_fits": fits,
            "V1_worst_rel": worst,
            "S_f_max_abs": mp.nstr(max(abs(S_f[n]) for n in S_f), 12),
            "S_f_max_abs_n": int(max(S_f, key=lambda k: abs(S_f[k]))),
        },
        "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())