"""nbval.py -- validation of the Vasyunin autocorrelation formula (eq:vasyunin) and of the
linear functional c_n = (gamma-1-log n)/n against direct numerical integration on the
critical line (Euler-Maclaurin zeta, truncation T; agreement = size of the 1/T tail)."""
import numpy as np
from math import gcd, log, pi
GAMMA = 0.5772156649015328606
def _zh_chunk(t):
    s = 0.5 + 1j*t
    M = int(np.max(np.abs(t))) + 60
    n = np.arange(1, M); logn = np.log(n)
    z = np.sum(np.exp(-s[:,None]*logn[None,:]), axis=1)
    z += M**(1-s)/(s-1) + 0.5*M**(-s)
    z += s*M**(-s-1.0)/12 - s*(s+1)*(s+2)*M**(-s-3.0)/720
    return z
def zeta_half(t):
    t = np.asarray(t, dtype=float); out = np.empty(len(t), dtype=complex)
    for i in range(0, len(t), 4000): out[i:i+4000] = _zh_chunk(t[i:i+4000])
    return out
def Smix(m, n):
    if n == 1: return 0.0
    k = np.arange(1, n)
    return float(np.sum(((m*k) % n)/n / np.tan(pi*k/n)))
def A(m, n):
    g = gcd(m, n); m, n = m//g, n//g
    out = (log(2*pi)-GAMMA)/2*(m+n)
    if m != n: out += (m-n)/2*log(n/m)
    out -= pi/2*(Smix(m, n) + Smix(n, m))
    return out/(m*n)/g
if __name__ == "__main__":
    T=6000; dt=0.002
    t=np.arange(dt/2, T, dt)
    z=zeta_half(t); z2=np.abs(z)**2
    w=z2/(0.25+t**2)*dt/(2*pi)*2
    print("pair  direct  formula   (tail ~ 1e-3)")
    for (m,n) in [(1,1),(1,2),(2,3),(1,3),(3,4),(5,6)]:
        direct = float(np.sum(w*np.cos(t*log(n/m))))
        print(f"({m},{n}): {direct:.6f}  {np.sqrt(m*n)*A(m,n):.6f}")
    for n in [1,2,3]:
        cd = 2*float(np.sum(np.real(z*np.exp(-1j*t*log(n)))/(0.25+t**2)))*dt/(2*pi)/np.sqrt(n)
        print(f"c_{n}: direct {cd:.6f}  formula {(GAMMA-1-log(n))/n:.6f}")