E-009-suzuki-screw-function.py

{
  "experiment": "E-009a",
  "task": "3.4c (D-018 step 1): Suzuki screw function g(t)",
  "description": "mpmath implementation of Suzuki's explicit screw function g(t) (arXiv:2606.09096v1 eq. 1.3); validated against eq. 2.2 near-origin expansion; tabulated on [0,10] step 1e-3 for E-009b Galerkin",
  "script": "research/numerics/scripts/E-009-suzuki-screw-function.py",
  "parameters": {
    "t_max": 10.0,
    "step": 0.001,
    "validation_dps": 30,
    "n_sieve_max": 22028
  },
  "checks": {
    "A1_lerch_convention": true,
    "A2_evenness": true,
    "A3_A_limit": true,
    "A4_r_curvature_converges": true,
    "B_table_finite": true,
    "B_spot_check_vs_mpmath": true
  },
  "key_values": {
    "A_eq2.2": 0.707546365655439,
    "c2_lim_r_over_t2": -0.8750599754042412,
    "digamma_coeff_-half(psi(1/4)-log pi)": 2.6860917096128327,
    "Phi(1,2,1/4)": 17.19732915450711,
    "g(log2)": -0.06355558603756917,
    "g(1)": -0.0440073052368517,
    "g(5)": -0.045371349558830154,
    "g(10)": -0.06209008245757719
  },
  "caveats": [
    "float64 table: Lerch combination has cancellation for small t (~4 digits lost at t=1e-3); spot-checked vs mpmath to 1e-9",
    "prime sum is piecewise linear in t (kinks at log n) — exact via cumsum, no interpolation error",
    "g even: only t >= 0 tabulated"
  ],
  "elapsed_seconds": 3.0,
  "timestamp_utc": "2026-08-18T01:44:26.674491+00:00",
  "row_count": 10001
}

"""E-009a — Task 3.4c (Q-3.4c, D-018 plan step 1): Suzuki's screw function g(t).

Suzuki (arXiv:2606.09096v1, eq. 1.3) — explicit screw function of zeta:

  g(t) = -4(e^{t/2} + e^{-t/2} - 2)
         + Sum_{n <= e^{|t|}} (Lambda(n)/sqrt(n)) (|t| - log n)
         - (|t|/2)(psi(1/4) - log pi)
         - (1/4)[ Phi(1,2,1/4) - e^{-|t|/2} Phi(e^{-2|t|}, 2, 1/4) ]

(psi = digamma, Phi = Hurwitz-Lerch transcendent Phi(z,s,a) =
Sum_{k>=0} z^k/(k+a)^s; Lambda = von Mangoldt; prime sum = 0 for
|t| < (1/2) log 2.)

Near-origin expansion (Suzuki eq. 2.2):
  g(t) = (1/2)|t| log|t| + A|t| + (prime sum) + r(t),
  r even C^2, r = O(t^2),
  A = (1/2)(log(2pi) - psi(2)) = (1/2)(log(2pi) + gamma_E - 1)
    = 0.7075463656...

g is a Krein-Langer screw function IFF RH (this is the Cor 1.6 bridge).

Part A (validation, mpmath dps 30):
  A1. Lerch convention: Phi(1,2,1/4) == psi'(1/4) = polygamma(1, 1/4).
  A2. Evenness: g(-t) == g(t) at t = 0.7, 2.3, 5.1 (exercises prime sum).
  A3. A-limit: g(t)/t - (1/2) log t -> A as t -> 0+ (prime sum = 0 there);
      residual must equal r(t)/t = c2*t + o(t) with c2 from A4.
  A4. r-curvature: r(t)/t^2 -> c2 (constant) as t -> 0+, t <= 0.25.

Part B (tabulation, float64 + mpmath hybrid):
  g(t) on t in [0, 10], step 1e-3 (mesh for the E-009b Galerkin stage;
  g is even so t >= 0 suffices). Prime-sum part via exact cumsum arrays
  (P(t) piecewise linear in t with kinks at log n). Lerch combination:
  direct numpy series for t > 0.35 (z = e^{-2t} < 0.5, fast), mpmath
  lerchphi for t <= 0.35 (slow-converging near z = 1).

Outputs: research/numerics/data/E-009-suzuki-screw-function.csv (+ .json)
"""
import json
import sys
import time
from datetime import datetime, timezone
from pathlib import Path

import numpy as np
import mpmath as mp

REPO = Path(__file__).resolve().parents[3]
DATA = REPO / "research" / "numerics" / "data"
CSV = DATA / "E-009-suzuki-screw-function.csv"
JSON = DATA / "E-009-suzuki-screw-function.json"

T_MAX = 10.0          # tabulation range [0, T_MAX] (g even)
H = 1e-3              # mesh step
DPS = 30              # validation precision
A_CONST = 0.5 * (np.log(2 * np.pi) + np.euler_gamma - 1.0)  # eq. 2.2 constant


# ---------------------------------------------------------------- von Mangoldt
def mangoldt_sieve(n_max: int) -> np.ndarray:
    """lam_log[n] = log p if n = p^k (prime power), else 0."""
    lam_log = np.zeros(n_max + 1)
    is_comp = np.zeros(n_max + 1, dtype=bool)
    for p in range(2, n_max + 1):
        if not is_comp[p]:
            if p * p <= n_max:
                is_comp[p * p :: p] = True
            pk = p
            while pk <= n_max:
                lam_log[pk] = np.log(p)
                pk *= p
    return lam_log


N_MAX = int(np.exp(T_MAX)) + 2
LAM_LOG = mangoldt_sieve(N_MAX)
NS = np.arange(2, N_MAX + 1)
LOGS = np.log(NS)
LAM = LAM_LOG[2:] / np.sqrt(NS)          # Lambda(n)/sqrt(n)
C1 = np.cumsum(LAM)                       # Sum_{n<=x} Lambda(n)/sqrt(n)
C0 = np.cumsum(LAM * LOGS)                # Sum_{n<=x} Lambda(n) log(n)/sqrt(n)


def prime_sum_vec(t: np.ndarray) -> np.ndarray:
    """P(t) = Sum_{n <= e^t} (Lambda(n)/sqrt(n))(t - log n), vectorized."""
    m = np.searchsorted(LOGS, t, side="right")     # count of n with log n <= t
    out = np.zeros_like(t)
    pos = m > 0
    idx = m[pos] - 1
    out[pos] = C1[idx] * t[pos] - C0[idx]
    return out


def prime_sum_mp(t: mp.mpf) -> mp.mpf:
    s = mp.mpf(0)
    n = 2
    while mp.log(n) <= t:
        s += (LAM_LOG[n] / mp.sqrt(n)) * (t - mp.log(n))
        n += 1
    return s


# ------------------------------------------------------------------- g(t), mp
def g_mp(t) -> mp.mpf:
    t = mp.mpf(t)
    a = abs(t)
    s = -4 * (mp.exp(a / 2) + mp.exp(-a / 2) - 2)
    s += prime_sum_mp(a)
    s -= (a / 2) * (mp.digamma(mp.mpf(1) / 4) - mp.log(mp.pi))
    s -= mp.mpf(1) / 4 * (mp.lerchphi(1, 2, mp.mpf(1) / 4)
                          - mp.exp(-a / 2) * mp.lerchphi(mp.exp(-2 * a), 2,
                                                         mp.mpf(1) / 4))
    return s


# --------------------------------------------------------------- g(t), float
KK = np.arange(0, 80)
DENOM = (KK + 0.25) ** 2


def lerch_phi2_fast(z: np.ndarray) -> np.ndarray:
    """Phi(z,2,1/4) by direct series; valid for |z| <= 0.5 (t >= 0.35)."""
    powers = np.power.outer(z, KK)          # (n_z, K)
    return (powers / DENOM).sum(axis=1)


def g_float_table(t: np.ndarray) -> np.ndarray:
    g = -4.0 * (np.exp(t / 2) + np.exp(-t / 2) - 2.0)
    g += prime_sum_vec(t)
    dig_coeff = -0.5 * float(mp.digamma(mp.mpf(1) / 4) - mp.log(mp.pi))
    g += dig_coeff * t
    phi1 = float(mp.lerchphi(1, 2, mp.mpf(1) / 4))   # = polygamma(1, 1/4)
    lerch_comb = np.empty_like(t)
    slow = t <= 0.35
    if slow.any():
        mp.dps = 25
        for i in np.flatnonzero(slow):
            ti = mp.mpf(float(t[i]))
            val = (mp.lerchphi(1, 2, mp.mpf(1) / 4)
                   - mp.exp(-ti / 2) * mp.lerchphi(mp.exp(-2 * ti), 2,
                                                   mp.mpf(1) / 4))
            lerch_comb[i] = float(val) / -4.0
        mp.dps = DPS
    if (~slow).any():
        z = np.exp(-2.0 * t[~slow])
        lerch_comb[~slow] = -0.25 * (phi1 - np.exp(-t[~slow] / 2)
                                     * lerch_phi2_fast(z))
    g += lerch_comb
    g[t == 0.0] = 0.0
    return g


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

    # ---- A1: Lerch convention
    phi1_mp = mp.lerchphi(1, 2, mp.mpf(1) / 4)
    trig1_mp = mp.polygamma(1, mp.mpf(1) / 4)
    a1_err = abs(phi1_mp - trig1_mp)
    check_a1 = bool(a1_err < mp.mpf("1e-25"))

    # ---- A2: evenness (t values exercising the prime sum)
    ev_ts = [mp.mpf("0.7"), mp.mpf("2.3"), mp.mpf("5.1")]
    ev_err = max(abs(g_mp(t) - g_mp(-t)) for t in ev_ts)
    check_a2 = bool(ev_err < mp.mpf("1e-20"))

    # ---- A4: r(t)/t^2 -> c2 (t <= 0.25: prime sum = 0)
    A_mp = (mp.log(2 * mp.pi) + mp.euler - 1) / 2
    r_ts = [mp.mpf("0.25"), mp.mpf("0.125"), mp.mpf("0.0625"),
            mp.mpf("0.03125"), mp.mpf("0.015625")]
    r_ratios = []
    for t in r_ts:
        r = g_mp(t) - (t / 2) * mp.log(t) - A_mp * t     # prime sum = 0 here
        r_ratios.append(r / t**2)
    c2 = r_ratios[-1]
    r_gap_last = abs(r_ratios[-1] - r_ratios[-2])
    r_gap_prev = abs(r_ratios[-2] - r_ratios[-3])
    check_a4 = bool(r_gap_last < r_gap_prev)  # gaps shrinking => converging

    # ---- A3: g(t)/t - (1/2) log t -> A; residual ~ c2 * t
    a3_ts = [mp.mpf("1e-2"), mp.mpf("1e-3"), mp.mpf("1e-4")]
    a3_resid = [g_mp(t) / t - mp.log(t) / 2 - A_mp for t in a3_ts]
    # expected residual ~ c2*t; test relative agreement within 25%
    a3_rel = [abs((res - c2 * t) / (c2 * t)) if c2 != 0 else abs(res)
              for res, t in zip(a3_resid, a3_ts)]
    check_a3 = bool(a3_rel[-1] < 0.25 and abs(a3_resid[-1]) < mp.mpf("3e-3"))

    print(f"A1 lerch==polygamma(1,1/4): err {mp.nstr(a1_err, 5)}")
    print(f"A2 evenness max err: {mp.nstr(ev_err, 5)}")
    print(f"A4 r(t)/t^2 -> c2 = {mp.nstr(c2, 12)} "
          f"(gaps {mp.nstr(r_gap_prev, 3)} -> {mp.nstr(r_gap_last, 3)})")
    print(f"A3 residuals g(t)/t - (1/2)log t - A at 1e-2/1e-3/1e-4: "
          f"{[mp.nstr(x, 6) for x in a3_resid]}  (c2*t predicted)")

    # ---- Part B: table on [0, T_MAX]
    t_mesh = np.arange(0.0, T_MAX + 0.5 * H, H)
    g_vals = g_float_table(t_mesh)
    p_vals = prime_sum_vec(t_mesh)
    check_b = bool(np.all(np.isfinite(g_vals)))
    # spot check: float table vs mpmath at a few points
    spot_ts = [0.123, 1.0, 3.7, 7.0, 10.0]
    spot_err = max(abs(g_mp(x) - mp.mpf(float(g_vals[int(round(x / H))])))
                   for x in spot_ts)
    check_spot = bool(spot_err < mp.mpf("1e-9"))
    print(f"Part B: {len(t_mesh)} points, spot-check err vs mpmath: "
          f"{mp.nstr(spot_err, 3)}")

    checks = {
        "A1_lerch_convention": check_a1,
        "A2_evenness": check_a2,
        "A3_A_limit": check_a3,
        "A4_r_curvature_converges": check_a4,
        "B_table_finite": check_b,
        "B_spot_check_vs_mpmath": check_spot,
    }

    DATA.mkdir(parents=True, exist_ok=True)
    with CSV.open("w", encoding="utf-8", newline="") as f:
        f.write("t,g_t,prime_sum_t\n")
        for ti, gi, pi in zip(t_mesh, g_vals, p_vals):
            f.write(f"{ti:.4f},{gi:.12e},{pi:.12e}\n")

    meta = {
        "experiment": "E-009a",
        "task": "3.4c (D-018 step 1): Suzuki screw function g(t)",
        "description": "mpmath implementation of Suzuki's explicit screw "
                       "function g(t) (arXiv:2606.09096v1 eq. 1.3); "
                       "validated against eq. 2.2 near-origin expansion; "
                       "tabulated on [0,10] step 1e-3 for E-009b Galerkin",
        "script": "research/numerics/scripts/E-009-suzuki-screw-function.py",
        "parameters": {"t_max": T_MAX, "step": H, "validation_dps": DPS,
                       "n_sieve_max": N_MAX},
        "checks": checks,
        "key_values": {
            "A_eq2.2": float(A_CONST),
            "c2_lim_r_over_t2": float(c2),
            "digamma_coeff_-half(psi(1/4)-log pi)":
                float(-(mp.digamma(mp.mpf(1) / 4) - mp.log(mp.pi)) / 2),
            "Phi(1,2,1/4)": float(phi1_mp),
            "g(log2)": float(g_mp(mp.log(2))),
            "g(1)": float(g_mp(1)),
            "g(5)": float(g_mp(5)),
            "g(10)": float(g_vals[-1]),
        },
        "caveats": [
            "float64 table: Lerch combination has cancellation for small t "
            "(~4 digits lost at t=1e-3); spot-checked vs mpmath to 1e-9",
            "prime sum is piecewise linear in t (kinks at log n) — exact "
            "via cumsum, no interpolation error",
            "g even: only t >= 0 tabulated",
        ],
        "elapsed_seconds": round(time.time() - t0, 1),
        "timestamp_utc": datetime.now(timezone.utc).isoformat(),
        "row_count": int(len(t_mesh)),
    }
    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())