E-021-diag3.py

"""Diag 3 — decisive: lambda_200 three independent ways + E-011 eta vs
diag eta at j = 90..200.
(1) zero side: lambda_200 = sum_rho [1-(1-1/rho)^200] (zeros cache)
(2) E-011 machinery eta (flint stieltjes + series division) -> S_f(200)
(3) diag eta (mpmath tau, dps 1700) -> S_f(200)
Truth: S_f(200) = S_inf(200) + 1 - lambda_200.
"""
import numpy as np
import mpmath as mp
import flint

zs = np.load("research/numerics/data/E-021-zeros-cache.npy")

# ---- (1) zero-side lambda_200
mp.mp.dps = 60
n = 200
lam = mp.mpf(0)
for g in zs:
    r = mp.mpc(0.5, float(g))
    lam += 2 * mp.re(1 - (1 - 1 / r) ** n)   # conjugate pair
# tail beyond T=5000: |1-(1-1/r)^n| <= 2n/|rho| for |rho| >> n (rough);
# more precisely contribution ~ n * sum 1/|rho|^2 tail ~ n * logT/(2pi T)
tail_est = n * np.log(5000 / (2 * np.pi)) / (2 * np.pi * 5000)
print(f"lambda_200 zero-side = {mp.nstr(lam, 15)}  (tail ~ {tail_est:.2e})")
print(f"  E-015 value        = 306.6557648513447664127226")

S_INF_200 = mp.mpf("303.515388171")   # E-021 CSV (12 digits shown; fine)

# ---- (2) E-011 eta via flint stieltjes + series division
flint.ctx.prec = 800
KMAX = 210
gam = [mp.mpf(flint.acb.stieltjes(k).real.str(80, radius=False))
       for k in range(KMAX + 2)]
mp.mp.dps = 150
N_ = [mp.mpf(0)] * (KMAX + 2)
D_ = [mp.mpf(0)] * (KMAX + 2)
N_[0] = D_[0] = 1
for m in range(0, KMAX + 1):
    D_[m + 1] = (-1) ** m * gam[m] / mp.factorial(m)
for m in range(1, KMAX + 1):
    N_[m + 1] = -(-1) ** m * gam[m] / mp.factorial(m - 1)
E = [mp.mpf(0)] * (KMAX + 2)
for k in range(0, KMAX + 2):
    s = N_[k]
    for i in range(1, k + 1):
        s -= D_[i] * E[k - i]
    E[k] = s
eta_e011 = [E[k + 1] for k in range(KMAX + 1)]

# ---- (3) diag eta (mpmath hurwitz tau, dps 1700)
rho = [mp.mpc(0.5, float(g)) for g in zs]
mp.mp.dps = 100
sigma = {1: sum(1 / (r * (1 - r)) for r in rho)}
for m in range(2, 17):
    sigma[m] = 2 * mp.re(sum(r ** (-m) for r in rho))
eta_diag = []
for j in range(0, KMAX + 1):
    mp.mp.dps = 1700
    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 = 100
    sig = sigma[j + 1] if j + 1 <= 16 else mp.mpf(0)
    eta_diag.append(tau + (-1) ** (j + 1) * sig + (-1) ** j)

print(f"\n{'j':>4} {'eta_e011 (flint)':>22} {'eta_diag (mpmath)':>22} {'diff':>12}")
for j in list(range(88, 100)) + list(range(190, 201)):
    mp.mp.dps = 100
    d = abs(eta_e011[j] - eta_diag[j])
    print(f"{j:>4} {mp.nstr(eta_e011[j], 12):>22} "
          f"{mp.nstr(eta_diag[j], 12):>22} {mp.nstr(d, 3):>12}")

# ---- S_f(200) both etas
for name, eta in (("e011", eta_e011), ("diag", eta_diag)):
    mp.mp.dps = 300
    s = mp.mpf(0)
    c = mp.mpf(1)
    for j in range(1, 201):
        c = c * (200 - j + 1) / j
        s += c * eta[j - 1]
    mp.mp.dps = 100
    print(f"S_f(200) [{name} eta] = {mp.nstr(s, 15)}")
print(f"S_f(200) truth (S_inf+1-lam, zero-side lam) = "
      f"{mp.nstr(S_INF_200 + 1 - lam, 12)}")
print(f"S_f(200) truth (S_inf+1-lam, E-015 lam)     = "
      f"{mp.nstr(S_INF_200 + 1 - mp.mpf('306.6557648513447664127226'), 12)}")