E-021-diag2.py

"""Diag 2: S_f(n) with three eta variants, find explosion law.
(i)  tau_j at dps 50 (current E-021 behavior)
(ii) tau_j at dps 1700
(iii) eta_j = (-1)^{j+1}[(1-2^{-(j+1)})zeta(j+1) - 1 + sigma_{j+1}]
      computed at dps 1700 (cancellation-tolerant: loses 0.48j digits
      of 1700 -> good to j ~ 3000)
"""
import numpy as np
import mpmath as mp

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

mp.mp.dps = 100
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 make_eta(kind):
    out = []
    for j in range(0, JMAX + 1):
        sig = sigma[j + 1] if j + 1 <= 16 else mp.mpf(0)
        if kind == "i":
            mp.mp.dps = 50
            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
            out.append(tau + (-1) ** (j + 1) * sig + (-1) ** j)
        elif kind == "ii":
            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
            out.append(tau + (-1) ** (j + 1) * sig + (-1) ** j)
        else:
            mp.mp.dps = 1700
            if j == 0:
                eta = (mp.mpf(-0.5) * (mp.log(4 * mp.pi) + mp.euler)
                       - sig + 1)
            else:
                small = ((1 - mp.mpf(2) ** (-(j + 1))) * mp.zeta(j + 1) - 1)
                eta = (-1) ** (j + 1) * (small + sig)
            mp.mp.dps = 100
            out.append(eta)
    return out

def sf_curve(eta, ns, dps):
    res = {}
    for n in ns:
        mp.mp.dps = dps
        s = mp.mpf(0)
        c = mp.mpf(1)
        for j in range(1, n + 1):
            c = c * (n - j + 1) / j
            s += c * eta[j - 1]
        res[n] = +s
    mp.mp.dps = 100
    return res

NS = [100, 200, 300, 400, 500, 700, 1000]
eta_i = make_eta("i")
eta_ii = make_eta("ii")
eta_iii = make_eta("iii")

sf_i = sf_curve(eta_i, NS, 500)
sf_ii = sf_curve(eta_ii, NS, 500)
sf_iii = sf_curve(eta_iii, NS, 500)

print(f"{'n':>5} {'S_f (i) dps50':>22} {'S_f (ii) dps1700':>22} "
      f"{'S_f (iii) direct':>22}")
for n in NS:
    print(f"{n:>5} {mp.nstr(sf_i[n], 12):>22} {mp.nstr(sf_ii[n], 12):>22} "
          f"{mp.nstr(sf_iii[n], 12):>22}")

# eta differences
print("\neta variant diffs (max |i-iii|, max |ii-iii|) over j ranges:")
for lo, hi in ((0, 16), (16, 100), (100, 300), (300, 700), (700, 1201)):
    d1 = max(abs(eta_i[j] - eta_iii[j]) for j in range(lo, hi))
    d2 = max(abs(eta_ii[j] - eta_iii[j]) for j in range(lo, hi))
    big = max(range(lo, hi), key=lambda j: abs(eta_i[j] - eta_iii[j]))
    print(f"  j {lo:>4}-{hi:>4}: i-iii {mp.nstr(d1, 3):>10} "
          f"(worst j={big})   ii-iii {mp.nstr(d2, 3):>10}")