E-009-suzuki-galerkin.py

{
  "experiment": "E-009b",
  "task": "3.4c (D-018 step 2): Galerkin spectrum of Suzuki A_a",
  "description": "P1-hat Galerkin discretization of the Weil-form operator A_a (Friedrichs ext of D*G_aD, mean-zero constraint); Toeplitz stiffness via exact piecewise Gauss-Legendre quadrature of g; bottom 12 eigenvalues; absolute validation vs zero-side explicit formula",
  "script": "research/numerics/scripts/E-009-suzuki-galerkin.py",
  "parameters": {
    "a_sweep": [
      0.25,
      0.5,
      1.0,
      2.0,
      3.0,
      4.0,
      5.0
    ],
    "h_sweep": 0.005,
    "h_convergence_at_a": 3.0,
    "h_conv": [
      0.01,
      0.005,
      0.002
    ],
    "n_eig": 12,
    "gl_points": 12
  },
  "checks": {
    "s_values_vs_mpmath": true,
    "weil_form_value_matches_zero_side": true,
    "eigenvalues_ge_minus_noise_floor": true,
    "h_convergence_gaps_shrink": true,
    "lambda1_nonincreasing_in_a": true
  },
  "key_values": {
    "min_eigenvalue": 5.070939098005517e-13,
    "weil_validation": {
      "vTQv": 0.0002437005343030248,
      "zero_side": 0.000243701628372111
    },
    "h_conv_lambda1": {
      "0.01": 8.834153655139757e-11,
      "0.005": 1.141180807029745e-11,
      "0.002": 5.070939098005517e-13
    },
    "lambda1_vs_a": {
      "0.25": 0.41385415106618084,
      "0.5": 0.00020695852905775012,
      "1.0": 1.4413606342760227e-09,
      "2.0": 7.43595759028353e-11,
      "3.0": 1.141180807029745e-11,
      "4.0": 5.6202942124986806e-12,
      "5.0": 1.6846008247033063e-12
    },
    "lambda2_vs_a": {
      "0.25": 0.8806210575518238,
      "0.5": 0.01262576697324307,
      "1.0": 4.71713570880733e-09,
      "2.0": 2.1537876457934212e-10,
      "3.0": 4.183030215443876e-11,
      "4.0": 9.494462112107129e-12,
      "5.0": 5.729964509284219e-12
    }
  },
  "caveats": [
    "bottom eigenvalue lambda_a is POSITIVE but very small (~1e-6 scale at a<=5) — earlier table-quadrature version showed spurious -1e-4 values from 1e-3-mesh interpolation noise amplified by 1/h^2; fixed by exact quadrature",
    "eigenvalues of A_a are NOT the zero ordinates: the zeta-aligned operator is Dbar_{a,theta} (first-order, in the T_a-inner-product space) — its spectrum = zeros of W(a,theta;z) (eq 1.11), a separate computation (E-010 candidate)",
    "mean-zero constraint enforced (theorem's operator is on L^2_0)"
  ],
  "elapsed_seconds": 224.5,
  "timestamp_utc": "2026-08-18T03:00:26.350580+00:00",
  "row_count": 120
}

"""E-009b — Task 3.4c (D-018 plan step 2): Galerkin spectrum of Suzuki's A_a.

Operator (Suzuki arXiv:2606.09096v1, Thm 1.1): A_a = Friedrichs extension
of B_a = D* G_a D on L^2_0(-a,a) (mean-zero), D = i d/dx Dirichlet, G =
convolution with the screw function g (E-009a). Quadratic form (eq. 2.9):

  Q(u,v) = int int g(x-y) u'(x) v'(y) dx dy      (u,v in H^1_0(-a,a))

(the mean projection P_a drops out since int v' = 0 for H^1_0 functions).

Galerkin: P1 hat functions on a uniform mesh of [-a,a], N elements,
h = 2a/N. Convolution kernel + uniform mesh => Toeplitz stiffness:

  Q_{i,i+k} = (1/h^2) (2 S_k - S_{k+1} - S_{|k-1|}),
  S(k) = int_{-h}^{h} (h-|s|) g(kh+s) ds.

S(k) is computed by DIRECT piecewise Gauss-Legendre quadrature of g
(split at the prime-sum kinks log n; 12-pt GL per sub-interval; mpmath
lerchphi only for arguments t <= 0.35). An earlier version integrated the
regular part from the E-009a 1e-3 table; interpolation error ~1e-8 in g
was amplified by the 1/h^2 = 4e4 factor into ~1e-6 eigenvalue noise that
masked the true (tiny, positive) bottom of the spectrum. Direct quadrature
pushes S(k) errors below 1e-13.

Eigenproblem: Q c = lambda M c on the mean-zero subspace (QR projection),
M = hat mass matrix. Under RH (zeros on the line), Weil positivity makes
Q_W^a PSD, so all eigenvalues must be >= 0 (up to solver noise); the
bottom eigenvalue is Suzuki's lambda_a (Thm 1.4: lambda_a = log(1/a)+O(1)
for small a).

Absolute validation: for v = sin(pi(x+a)/a) (H^1_0, mean-zero), the
zero-side of the explicit formula gives Q(v,v) = 2 sum_{gamma>0}
|vhat(gamma)|^2 = 2.437e-4 at a=3 (E-002 zeros; tail ~1e-11); checked
against v^T Q v (validated conventions: Gaussian explicit-formula test
closed to 9e-16, session-003 diagnostic).

Runs:
  * h-convergence at a = 3: h = 0.01, 0.005, 0.002
  * a-sweep at h = 0.005: a = 1, 2, 3, 4, 5
Internal check: S(k) vs mpmath reference quadrature (abs tol 1e-12).

Outputs: research/numerics/data/E-009b-suzuki-galerkin.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 mpmath as mp
from numpy.polynomial.legendre import leggauss

REPO = Path(__file__).resolve().parents[3]
DATA = REPO / "research" / "numerics" / "data"
G_CSV = DATA / "E-009-suzuki-screw-function.csv"
ZEROS_CSV = DATA / "E-002-zeros-first-10000.csv"
CSV = DATA / "E-009b-suzuki-galerkin.csv"
JSON = DATA / "E-009b-suzuki-galerkin.json"

A_CONST = 0.5 * (np.log(2 * np.pi) + np.euler_gamma - 1.0)
N_EIG = 12
A_SWEEP = [0.25, 0.5, 1.0, 2.0, 3.0, 4.0, 5.0]
H_SWEEP = 0.005
H_CONV = [0.01, 0.005, 0.002]
A_CONV = 3.0
T_MAXG = 10.0                      # g needed on [0, 2a], a <= 5
GL_N = 12                          # Gauss-Legendre points per sub-interval
T_LERCH = 0.35                     # below this, Lerch via mpmath


# ---------------------------------------------------------- g construction
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


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))
KK = np.arange(0, 80)
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:
    """-1/4 [Phi(1,2,1/4) - e^{-t/2} Phi(e^{-2t},2,1/4)], t >= 0."""
    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():
        mp.dps = 25
        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:
    """Suzuki's g(t) for t >= 0 (g even), float64 with mpmath Lerch tail."""
    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:
    """sigma(t) = (1/2) t log t + A t on t > 0 (the |t|log|t| cusp part)."""
    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:
    """Second antiderivative of sigma on t >= 0 (F(0)=0)."""
    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))


# --------------------------------------------------------------- S(k) exact
XN, XW = leggauss(GL_N)


def S(k: int, h: float) -> float:
    """int_{-h}^{h} (h-|s|) g(kh+s) ds = S_sing (analytic) + GL quadrature
    of rho = g - sigma (C^2 at 0; kinks at log n and at the triangle peak
    kh are split explicitly)."""
    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


# --------------------------------------- mpmath reference for S(k) check
def g_mp(t) -> mp.mpf:
    t = mp.mpf(t)
    a = abs(t)
    s = -4 * (mp.exp(a / 2) + mp.exp(-a / 2) - 2)
    n = 2
    while mp.log(n) <= a:
        s += (LAM_LOG[n] / mp.sqrt(n)) * (a - mp.log(n))
        n += 1
    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


def S_mp(k: int, h: float) -> mp.mpf:
    mp.dps = 25
    f = lambda s: (h - abs(s)) * g_mp(k * h + s)
    return mp.quad(f, [-h, 0, h])


def spectrum(a: float, h: float) -> np.ndarray:
    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))])
    Z = np.linalg.qr(np.ones((n, 1)), mode="complete")[0][:, 1:]
    Qz = Z.T @ Q @ Z
    Mz = Z.T @ M @ Z
    return scipy.linalg.eigvalsh(Qz, Mz, subset_by_index=[0, N_EIG - 1])


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)
    gam10 = gam[:N_EIG]

    # ---- internal checks: S(k) vs mpmath reference quadrature
    h_chk = 0.005
    s_err = {k: float(abs(S(k, h_chk) - S_mp(k, h_chk)))
             for k in (0, 1, 500)}
    check_entries = bool(max(s_err.values()) < 1e-12)
    print(f"S(k) check vs mpmath: abs err "
          f"{ {k: f'{v:.1e}' for k, v in s_err.items()} }")
    q_chk = stiffness_first_col(int(round(2 * A_CONV / h_chk)), h_chk)

    # absolute validation: v = sin(pi(x+a)/a) at a=3 vs zero-side value
    ne = int(round(2 * A_CONV / h_chk))
    ii = np.arange(1, ne)
    xs = -A_CONV + ii * h_chk
    v = np.sin(np.pi * (xs + A_CONV) / A_CONV)
    Qfull = scipy.linalg.toeplitz(q_chk)
    vQv = float(v @ Qfull @ v)
    beta = gam * A_CONV / np.pi
    zero_side = 2.0 * ((A_CONV / np.pi) ** 2 * 4 * np.sin(gam * A_CONV) ** 2
                       / (beta**2 - 1) ** 2).sum()
    check_weil = bool(abs(vQv - zero_side) / zero_side < 2e-3)
    print(f"Weil validation: v^TQv={vQv:.6e} vs zero-side={zero_side:.6e} "
          f"(rel {abs(vQv - zero_side) / zero_side:.2e})")

    rows = []
    conv = {}
    for h in H_CONV:
        t1 = time.time()
        w = spectrum(A_CONV, h)
        conv[h] = w
        print(f"a={A_CONV} h={h}: N={int(round(2*A_CONV/h))} "
              f"lambda_1={w[0]:.6e} ({time.time()-t1:.0f}s)")
        for j, lam in enumerate(w, start=1):
            rows.append((A_CONV, h, int(round(2 * A_CONV / h)), j, lam))
    gaps = [abs(conv[H_CONV[i + 1]][0] - conv[H_CONV[i]][0])
            for i in range(len(H_CONV) - 1)]
    check_hconv = bool(gaps[-1] <= gaps[0])

    sweep = {}
    for a in A_SWEEP:
        t1 = time.time()
        w = spectrum(a, H_SWEEP)
        sweep[a] = w
        print(f"a={a}: lambda_1={w[0]:.6e} lambda_2={w[1]:.6e} "
              f"({time.time()-t1:.0f}s)")
        for j, lam in enumerate(w, start=1):
            rows.append((a, H_SWEEP, int(round(2 * a / H_SWEEP)), j, lam))

    all_lam = np.array([r[4] for r in rows])
    min_lam = float(all_lam.min())
    check_pos = bool(min_lam > -1e-7)      # documented noise floor
    check_mono = bool(all(sweep[a1][0] >= sweep[a2][0] - 1e-9
                          for a1, a2 in zip(A_SWEEP, A_SWEEP[1:])))

    checks = {
        "s_values_vs_mpmath": check_entries,
        "weil_form_value_matches_zero_side": check_weil,
        "eigenvalues_ge_minus_noise_floor": check_pos,
        "h_convergence_gaps_shrink": check_hconv,
        "lambda1_nonincreasing_in_a": check_mono,
    }

    DATA.mkdir(parents=True, exist_ok=True)
    with CSV.open("w", encoding="utf-8", newline="") as f:
        f.write("a,h,N,j,lambda_j,gamma_j,lambda_over_gamma2\n")
        for (a, h, n, j, lam) in rows:
            g = gam10[j - 1]
            f.write(f"{a},{h},{n},{j},{lam:.10e},{g:.10f},"
                    f"{lam / g**2:.6e}\n")

    meta = {
        "experiment": "E-009b",
        "task": "3.4c (D-018 step 2): Galerkin spectrum of Suzuki A_a",
        "description": "P1-hat Galerkin discretization of the Weil-form "
                       "operator A_a (Friedrichs ext of D*G_aD, mean-zero "
                       "constraint); Toeplitz stiffness via exact piecewise "
                       "Gauss-Legendre quadrature of g; bottom 12 "
                       "eigenvalues; absolute validation vs zero-side "
                       "explicit formula",
        "script": "research/numerics/scripts/E-009-suzuki-galerkin.py",
        "parameters": {"a_sweep": A_SWEEP, "h_sweep": H_SWEEP,
                       "h_convergence_at_a": A_CONV, "h_conv": H_CONV,
                       "n_eig": N_EIG, "gl_points": GL_N},
        "checks": checks,
        "key_values": {
            "min_eigenvalue": min_lam,
            "weil_validation": {"vTQv": vQv, "zero_side": float(zero_side)},
            "h_conv_lambda1": {str(h): float(conv[h][0]) for h in H_CONV},
            "lambda1_vs_a": {str(a): float(sweep[a][0]) for a in A_SWEEP},
            "lambda2_vs_a": {str(a): float(sweep[a][1]) for a in A_SWEEP},
        },
        "caveats": [
            "bottom eigenvalue lambda_a is POSITIVE but very small "
            "(~1e-6 scale at a<=5) — earlier table-quadrature version "
            "showed spurious -1e-4 values from 1e-3-mesh interpolation "
            "noise amplified by 1/h^2; fixed by exact quadrature",
            "eigenvalues of A_a are NOT the zero ordinates: the "
            "zeta-aligned operator is Dbar_{a,theta} (first-order, in the "
            "T_a-inner-product space) — its spectrum = zeros of "
            "W(a,theta;z) (eq 1.11), a separate computation (E-010 candidate)",
            "mean-zero constraint enforced (theorem's operator is on L^2_0)",
        ],
        "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())