import numpy as np, json, time
from math import gcd, log, pi
GAMMA = 0.5772156649015328606
LOG2PI = log(2*pi)
SMEMO={}
def Smix(m, n):
    if n == 1: return 0.0
    key=(m % n, n)
    if key in SMEMO: return SMEMO[key]
    k = np.arange(1, n)
    v = float(np.sum(((m*k) % n)/n / np.tan(pi*k/n)))
    SMEMO[key]=v
    return v
C0 = (LOG2PI-GAMMA)/2
def Fexpr(r, b):
    out = C0*(b+r)
    if r != b: out += (b-r)/2*log(r/b)
    out -= pi/2*(Smix(b, r) + Smix(r, b))
    return out/(b*r)
resd=[]
for N in [5,10,20,40,80,160,320,500,640,800]:
    ns=list(range(1,N+1))
    Gf=np.zeros((N,N))
    for i,m in enumerate(ns):
        for j in range(i,N):
            n=ns[j]; g=gcd(m,n)
            Gf[i,j]=Gf[j,i]=Fexpr(n//g, m//g)/g
    cf=np.array([(GAMMA-1-log(n))/n for n in ns])
    w,V=np.linalg.eigh(Gf)
    thr=max(1e-13*w.max(),0.0)
    winv=np.where(w>thr,1/w,0.0)
    x=V@(winv*(V.T@cf))
    d2=1-float(cf@x)
    # reliability: residual and negative-eigenvalue mass
    res=float(np.linalg.norm(Gf@x-cf)/np.linalg.norm(cf))
    resd.append((N,d2,d2*log(N),float(w.min()),res))
    print(f"N={N:4d}: d2={d2:.6f}  d2*logN={d2*log(N):.5f}  min_eig={w.min():.2e}  rel_res={res:.1e}", flush=True)
json.dump(resd, open('nb_dN2.json','w'))
print("DONE", flush=True)