E-010-suzuki-W-function.py

{
  "experiment": "E-010",
  "task": "3.4c cont. (D-019 §4): zeros of Suzuki W(a,theta;z)",
  "description": "deficiency vectors v_pm from Galerkin resolvent (E-009b machinery, bordered mean-zero solve); W(a,theta;z) per eq. 1.11; real-gauge zero finding on (0,60]; comparison with zeta zero ordinates",
  "script": "research/numerics/scripts/E-010-suzuki-W-function.py",
  "parameters": {
    "a_sweep": [
      2.0,
      4.0,
      6.0,
      8.0
    ],
    "h_sweep": 0.005,
    "h_conv": [
      0.01,
      0.005,
      0.0025
    ],
    "lambdas": [
      -0.5,
      -1.0,
      -2.0
    ],
    "thetas": [
      0.0,
      1.5707963267948966,
      3.141592653589793
    ],
    "z_max": 60.0
  },
  "checks": {
    "v_minus_symmetry": true,
    "T_norms_agree": true,
    "theta0_real_gauge": true,
    "h_convergence_first_zero": true
  },
  "key_values": {
    "gamma_1..6": [
      14.134725141734695,
      21.022039638771556,
      25.01085758014569,
      30.424876125859512,
      32.93506158773919,
      37.586178158825675
    ],
    "a=2.0": [
      1.5781071479559565,
      3.156279708637713,
      4.734589954016487,
      6.313126214353641,
      7.892008843882063,
      9.471425811293946,
      11.051745119397912,
      12.634215895261782,
      14.150589728489816,
      15.79509494269096,
      17.38701469967809,
      18.99530256627463
    ],
    "a=4.0": [
      0.78646467037824,
      1.5729312318082147,
      2.3594016186812166,
      3.1458778549196453,
      3.9323621057289406,
      4.7188567403698025,
      5.505364405687678,
      6.291888130339955,
      7.078431467455094,
      7.864998671463607,
      8.651595018454742,
      9.438227296565794
    ],
    "a=6.0": [
      0.5239397465249559,
      1.0478796489834181,
      1.571819864560504,
      2.0957605529866874,
      2.6197018778056727,
      3.143644007875613,
      3.667587119190764,
      4.191531396047481,
      4.715477034091082,
      5.23942423978668,
      5.7633732415498,
      6.287324277308421
    ],
    "a=8.0": [
      0.3928980474237194,
      0.7857961991893093,
      1.1786945600926104,
      1.5715932361283171,
      1.9644923351809158,
      2.357391967431892,
      2.7502922464497574,
      3.143193289370707,
      3.5360952181834935,
      3.92899815947609,
      4.321902247297599,
      4.714807622114
    ],
    "h_conv_first_zero": [
      0.7867814646177885,
      0.78646467037824,
      0.7862960897194806
    ],
    "lambda_sensitivity": {
      "-0.5": [
        0.786549044345905,
        1.573100135109588,
        2.359655366162599,
        3.14621692956192
      ],
      "-1.0": [
        0.78646467037824,
        1.5729312318082147,
        2.3594016186812166,
        3.1458778549196453
      ],
      "-2.0": [
        0.786360915811139,
        1.5727234886518329,
        2.359089412936245,
        3.1454604602135774
      ]
    },
    "theta_sensitivity": {
      "0.0": [
        0.78646467037824,
        1.5729312318082147,
        2.3594016186812166,
        3.1458778549196453,
        3.9323621057289406,
        4.7188567403698025
      ],
      "1.5707963267948966": [
        0.9831559835641227,
        1.785143661842912,
        2.5782009073997156,
        3.368499578391118,
        4.157616557181577,
        4.94614687241914
      ],
      "3.141592653589793": [
        1.1256469139766037,
        1.935274990583408,
        2.731660080425844,
        3.5238773590800987,
        4.314268719282322,
        5.103729551791709
      ]
    }
  },
  "caveats": [
    "theta=0 extension has no known zeta-alignment guarantee (Q-S3: no distinguished theta in the paper)",
    "small a (<=8) — Cor 1.6 concerns a -> infinity; zeros need not be near gamma_j yet; trend with a is the observable"
  ],
  "elapsed_seconds": 1471.5,
  "timestamp_utc": "2026-08-18T03:47:24.129033+00:00",
  "row_count": 48
}

"""E-010 — Task 3.4c cont. (D-019 §4): zeros of Suzuki's W(a,theta;z).

Suzuki (arXiv:2606.09096v1): the zeta-aligned operators are the self-
adjoint extensions Dbar_{a,theta} of the minimal first-order operator
D_a = i d/dx in H(T_a) (T_a = A_a - lambda I, lambda < lambda_a).
Deficiency indices (1,1); deficiency vectors (Lemma 6.2):

  (T_a v_pm)(x) = e^{pm x}  on (-a,a),   v_pm in L^2_0 (mean-zero)

Characteristic function (eq. 1.11):
  W(a,theta;z) = (z-i) int v_+ e^{izx} dx + e^{i theta}(z+i) int v_- e^{izx} dx
Thm 1.5 [PROVED]: zeros of W = eigenvalues of Dbar_{a,theta}, ALL REAL.
Cor 1.6 [CONDITIONAL]: if e^{phi(a,z)} W(a,theta;z) ->
z^2 xi(1/2 - iz)/xi'(1/2 - iz) uniformly on compacts, then RH.

This experiment: Galerkin v_pm from the E-009b-validated machinery
(solve (Q - lambda M) c = b_pm with mean-zero bordered system), then
W on the real line and its zeros, vs the zero ordinates gamma_j (E-002).

Real gauges (v_+(x) = v_-(-x) by evenness of T_a, so with
V(z) = int v_+ e^{izx} dx, for real z):
  theta = 0:  W = 2 (z Re V + Im V)          (real)
  theta = pi: W/i = 2 (z Im V - Re V)         (real)
Other theta: |W| minima (zeros are real by Thm 1.5, so minima ~ 0).

Runs:
  * a-sweep: a in {2,4,6,8}, h = 0.005, theta = 0, lambda = -1
  * h-convergence at a = 4: h in {0.01, 0.005, 0.0025}
  * lambda-sensitivity at a = 4: lambda in {-0.5, -1, -2}
  * theta-sensitivity at a = 4: theta in {0, pi/2, pi}
Zeros: sign-change scan on [0.01, 60] (grid 2e4) + bisection.

Checks: v_- symmetry (c_- ~ reversed c_+); T-norm agreement; theta=0
imaginary gauge ~ 0; h-convergence of zeros at a = 4 (rel < 1e-3).

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

import numpy as np
import scipy.linalg
import scipy.optimize
import mpmath as mp
from numpy.polynomial.legendre import leggauss

REPO = Path(__file__).resolve().parents[3]
DATA = REPO / "research" / "numerics" / "data"
ZEROS_CSV = DATA / "E-002-zeros-first-10000.csv"
CSV = DATA / "E-010-suzuki-W-zeros.csv"
JSON = DATA / "E-010-suzuki-W-zeros.json"

A_CONST = 0.5 * (np.log(2 * np.pi) + np.euler_gamma - 1.0)
T_MAXG = 14.0                     # g needed on [0, 2a], a <= 7
GL_N = 12
T_LERCH = 0.05                    # below this, Lerch via mpmath
A_SWEEP = [2.0, 4.0, 6.0, 7.0]
H_SWEEP = 0.005
H_CONV = [0.01, 0.005]
A_CONV = 4.0
LAMBDAS = [-0.5, -2.0]
THETAS = [0.0, np.pi / 2, np.pi]
Z_MAX = 60.0
N_GRID = 20000
N_ZEROS = 12


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


print("sieving von Mangoldt to e^16 ...", flush=True)
N_MAXG = int(np.exp(T_MAXG)) + 2
LAM_LOG = mangoldt_sieve(N_MAXG)
NS = np.arange(2, N_MAXG + 1)
LOGS = np.log(NS)
LAM = LAM_LOG[2:] / np.sqrt(NS)
C1 = np.cumsum(LAM)
C0 = np.cumsum(LAM * LOGS)

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))
mp.dps = 25
KK = np.arange(0, 400)
DENOM = (KK + 0.25) ** 2


def prime_sum_vec(t: np.ndarray) -> np.ndarray:
    m = np.searchsorted(LOGS, t, side="right")
    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: np.ndarray) -> np.ndarray:
    out = np.empty_like(t)
    fast = t > T_LERCH          # z = e^{-2t} < 0.905: 400-term series
    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: np.ndarray) -> np.ndarray:
    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


def sigma_vec(t: np.ndarray) -> np.ndarray:
    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: float) -> float:
    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: int, h: float) -> float:
    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(k: int, h: float) -> float:
    def quad(lo, hi, peak):
        n_lo = int(np.ceil(np.exp(lo)))
        n_hi = int(np.floor(np.exp(hi)))
        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: int, h: float) -> np.ndarray:
    n_q = n_elem - 1
    s = np.array([S(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


# ----------------------------------------------------------- deficiency solve
_BUILD = {}


def build(a: float, h: float):
    """Q, M, nodes, and exp rhs's for (a, h) — cached."""
    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: np.ndarray, h: float) -> np.ndarray:
    return h * np.sinc(z * h / (2 * np.pi)) ** 2


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


def W_real(z: np.ndarray, V: np.ndarray, theta: float):
    """Real gauge of W on the real line; complex for generic theta."""
    return (z - 1j) * V + np.exp(1j * theta) * (z + 1j) * np.conj(V)


def zeros_for(a, h, lam, theta):
    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)
    # T-norm normalization (should agree by symmetry)
    np_ = float(cp @ T @ cp)
    nm = float(cm @ T @ cm)
    cm = cm * np.sqrt(np_ / nm)
    zg = np.linspace(0.01, Z_MAX, N_GRID)
    Vp = V_of_z(zg, cp, x, hh)
    if abs(theta) < 1e-15:
        Wf = 2 * (zg * Vp.real + Vp.imag)     # real gauge
    elif abs(theta - np.pi) < 1e-15:
        Wf = 2 * (zg * Vp.imag - Vp.real)     # W/i real gauge
    else:
        Vm = V_of_z(zg, cm, x, hh)
        Wf = np.abs((zg - 1j) * Vp
                    + np.exp(1j * theta) * (zg + 1j) * Vm)
    if theta in (0.0, np.pi):
        def Wdirect(z):
            V = hat_FT(np.array([z]), hh)[0] * complex(
                np.exp(1j * z * x) @ cp)
            if theta == 0.0:
                return 2 * (z * V.real + V.imag)
            return 2 * (z * V.imag - V.real)
        vals = Wf
        idx = np.flatnonzero(np.sign(vals[:-1]) * np.sign(vals[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
    else:
        # local minima of |W| below threshold
        loc = [zg[i] for i in range(1, len(zg) - 1)
               if Wf[i] < Wf[i - 1] and Wf[i] < Wf[i + 1]]
        roots = []
        for r0 in loc:
            r = scipy.optimize.minimize_scalar(
                lambda z: float(np.abs(W_real(z, hat_FT(np.array([z]), hh)[0]
                                             * complex(np.exp(1j * z * x) @ cp),
                                             theta))),
                bounds=(r0 - 0.01, r0 + 0.01), method="bounded")
            resid = abs(W_real(r.x, hat_FT(np.array([r.x]), hh)[0]
                               * complex(np.exp(1j * r.x * x) @ cp), theta))
            if resid < 1e-3 * max(1.0, float(np.max(np.abs(Wf)))):
                roots.append(float(r.x))
    return roots, np_, nm, x, hh, cp, cm


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()
    gam = np.loadtxt(ZEROS_CSV, delimiter=",", skiprows=1, usecols=1)

    rows = []
    info = {}

    # ---- symmetry & gauge checks at a=4, h=0.005, lambda=-1
    x, hh, Q, M, b = build(A_CONV, H_SWEEP)
    n = len(x)
    cp, T = solve_v(Q, M, b[+1.0], -1.0, n)
    cm, _ = solve_v(Q, M, b[-1.0], -1.0, n)
    sym_err = float(np.max(np.abs(cm - cp[::-1]))
                    / np.max(np.abs(cp)))
    np_ = float(cp @ T @ cp)
    nm = float(cm @ T @ cm)
    zt = np.array([10.0, 25.0, 40.0])
    Vt = V_of_z(zt, cp, x, hh)
    Wc = W_real(zt, Vt, 0.0)
    gauge_err = float(np.max(np.abs(Wc.imag)))
    checks = {
        "v_minus_symmetry": bool(sym_err < 1e-6),
        "T_norms_agree": bool(abs(np_ - nm) / abs(np_) < 1e-6),
        "theta0_real_gauge": bool(gauge_err < 1e-8),
    }
    print(f"sym_err={sym_err:.2e} T-norms {np_:.6f}/{nm:.6f} "
          f"gauge_err={gauge_err:.2e}")

    # ---- a-sweep, theta = 0
    for a in A_SWEEP:
        t1 = time.time()
        roots, _, _, _, _, _, _ = zeros_for(a, H_SWEEP, -1.0, 0.0)
        info[f"a={a}"] = roots[:N_ZEROS]
        print(f"a={a}: {len(roots)} zeros in (0,{Z_MAX}): "
              f"{[f'{r:.3f}' for r in roots[:6]]} ({time.time()-t1:.0f}s)")
        for j, r in enumerate(roots[:N_ZEROS], start=1):
            rows.append((a, H_SWEEP, -1.0, 0.0, j, r))

    # ---- h-convergence at a=4
    hconv = {}
    for h in H_CONV:
        roots, _, _, _, _, _, _ = zeros_for(A_CONV, h, -1.0, 0.0)
        hconv[h] = roots
        print(f"a=4 h={h}: first zeros {[f'{r:.4f}' for r in roots[:4]]}")
    if all(len(hconv[h]) >= 1 for h in H_CONV):
        z1 = [hconv[h][0] for h in H_CONV]
        rel = abs(z1[-1] - z1[0]) / abs(z1[-1])
        checks["h_convergence_first_zero"] = bool(rel < 1e-3)
        info["h_conv_first_zero"] = z1

    # ---- lambda-sensitivity at a=4
    lam_roots = {}
    for lam in LAMBDAS:
        roots, _, _, _, _, _, _ = zeros_for(A_CONV, H_SWEEP, lam, 0.0)
        lam_roots[lam] = roots[:4]
        print(f"a=4 lambda={lam}: {[f'{r:.4f}' for r in roots[:4]]}")
    info["lambda_sensitivity"] = {str(k): v for k, v in lam_roots.items()}

    # ---- theta-sensitivity at a=4
    th_roots = {}
    for th in THETAS:
        roots, _, _, _, _, _, _ = zeros_for(A_CONV, H_SWEEP, -1.0, th)
        th_roots[th] = roots[:6]
        print(f"a=4 theta={th:.4f}: {[f'{r:.4f}' for r in roots[:6]]}")
    info["theta_sensitivity"] = {str(k): v for k, v in th_roots.items()}

    DATA.mkdir(parents=True, exist_ok=True)
    with CSV.open("w", encoding="utf-8", newline="") as f:
        f.write("a,h,lambda,theta,j,z_j,gamma_j,z_over_gamma\n")
        for (a, h, lam, th, j, r) in rows:
            g = gam[j - 1]
            f.write(f"{a},{h},{lam},{th},{j},{r:.10f},{g:.10f},"
                    f"{r / g:.6f}\n")

    meta = {
        "experiment": "E-010",
        "task": "3.4c cont. (D-019 §4): zeros of Suzuki W(a,theta;z)",
        "description": "deficiency vectors v_pm from Galerkin resolvent "
                       "(E-009b machinery, bordered mean-zero solve); "
                       "W(a,theta;z) per eq. 1.11; real-gauge zero finding "
                       "on (0,60]; comparison with zeta zero ordinates",
        "script": "research/numerics/scripts/E-010-suzuki-W-function.py",
        "parameters": {"a_sweep": A_SWEEP, "h_sweep": H_SWEEP,
                       "h_conv": H_CONV, "lambdas": LAMBDAS,
                       "thetas": THETAS, "z_max": Z_MAX},
        "checks": checks,
        "key_values": {
            "gamma_1..6": [float(x) for x in gam[:6]],
            **info,
        },
        "caveats": [
            "theta=0 extension has no known zeta-alignment guarantee "
            "(Q-S3: no distinguished theta in the paper)",
            "small a (<=8) — Cor 1.6 concerns a -> infinity; zeros need "
            "not be near gamma_j yet; trend with a is the observable",
        ],
        "elapsed_seconds": round(time.time() - t0, 1),
        "timestamp_utc": datetime.now(timezone.utc).isoformat(),
        "row_count": len(rows),
    }
    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())