E-021-diag-eta.py

"""Diagnostic for E-021 V1 failure: compare eta_j conventions/precisions.
(a) E-021 zero-side, tau_j at dps 50 (current script behavior)
(b) E-021 zero-side, tau_j at dps 300
(c) E-011 Stieltjes series division (flint.acb.stieltjes) — validated
Then S_f(n) = sum C(n,j) eta_{j-1} (plain) with (c) at high dps for
n in {100, 500, 1000} vs E-021 CSV / E-015 lambda.
"""
import sys
import flint
import mpmath as mp
import numpy as np

JMAX = 60

# ---- (c) Stieltjes eta (E-011 machinery, serial)
mp.mp.dps = 120
flint.ctx.prec = 600
gam = []
for k in range(JMAX + 2):
    v = flint.acb(0, 0).stieltjes(k) if hasattr(flint.acb, "stieltjes") else None
    gam.append(v)
# flint stieltjes API: flint.acb.stieltjes? use arb via flint.arb? fallback:
try:
    g0 = flint.acb(1).stieltjes(0)
except Exception as e:
    print("stieltjes API probe failed:", e)
    sys.exit(1)

gam = [mp.mpf(flint.acb(1).stieltjes(k).real.str(30, radius=False))
       for k in range(JMAX + 2)]

def eta_from_stieltjes(gam, n_max):
    N_ = [mp.mpf(0)] * (n_max + 2)
    D_ = [mp.mpf(0)] * (n_max + 2)
    N_[0] = D_[0] = 1
    for m in range(0, n_max + 1):
        D_[m + 1] = (-1) ** m * gam[m] / mp.factorial(m)
    for m in range(1, n_max + 1):
        N_[m + 1] = -(-1) ** m * gam[m] / mp.factorial(m - 1)
    E = [mp.mpf(0)] * (n_max + 2)
    for k in range(0, n_max + 2):
        s = N_[k]
        for i in range(1, k + 1):
            s -= D_[i] * E[k - i]
        E[k] = s
    return [E[k + 1] for k in range(n_max + 1)]

eta_c = eta_from_stieltjes(gam, JMAX)

# ---- zeros + sigma sums (E-021 machinery)
zs = np.load("research/numerics/data/E-021-zeros-cache.npy")
mp.mp.dps = 80
rho = [mp.mpc(0.5, float(g)) for g in zs]
sigma = {}
for m in range(1, 17):
    if m == 1:
        sigma[1] = sum(1 / (r * (1 - r)) for r in rho)
    else:
        sigma[m] = 2 * mp.re(sum(r ** (-m) for r in rho))

def eta_zero(dps_tau):
    out = []
    for j in range(0, JMAX + 1):
        mp.mp.dps = dps_tau
        if j == 0:
            tau = mp.mpf(-0.5) * (mp.log(4 * mp.pi) + mp.euler)
        else:
            tau = mp.mpf(-0.5) ** (j + 1) * mp.zeta(j + 1, mp.mpf(0.5))
        mp.mp.dps = 80
        sig = sigma[j + 1] if j + 1 <= 16 else mp.mpf(0)
        out.append(tau + (-1) ** (j + 1) * sig + (-1) ** j)
    return out

eta_a = eta_zero(50)
eta_b = eta_zero(300)

print(f"{'j':>3} {'eta_c (stieltjes)':>22} {'eta_a (tau dps50)':>22} "
      f"{'eta_b (tau dps300)':>22}")
for j in range(0, JMAX + 1):
    mp.mp.dps = 80
    flag = ""
    if abs(eta_a[j] - eta_c[j]) > max(mp.mpf('1e-20'), 1e-6 * abs(eta_c[j])):
        flag = "  <-- a WRONG"
    if abs(eta_b[j] - eta_c[j]) > max(mp.mpf('1e-40'), 1e-6 * abs(eta_c[j])):
        flag += "  <-- b WRONG"
    print(f"{j:>3} {mp.nstr(eta_c[j], 12):>22} {mp.nstr(eta_a[j], 12):>22} "
          f"{mp.nstr(eta_b[j], 12):>22}{flag}")

# ---- S_f with eta_c at high dps
mp.mp.dps = 700
eta_big = [mp.mpf(mp.nstr(e, 60)) for e in eta_c]  # 60 digits is plenty here
# need eta up to n=1000 -> Stieltjes to k=1000 too expensive for diag;
# instead compare small n only where we trust: n <= 60
for n in (2, 10, 50, 60):
    s = mp.mpf(0)
    c = mp.mpf(1)
    for j in range(1, n + 1):
        c = c * (n - j + 1) / j
        s += c * eta_big[j - 1]
    print(f"S_f({n}) [plain, stieltjes eta] = {mp.nstr(s, 15)}")