import numpy as np, json, time
from math import gcd, log, pi

GAMMA = 0.5772156649015328606
LOG2PI = log(2*pi)
C_BD = 2 + GAMMA - log(4*pi)

SMEMO={}
def Smix(m, n):
    """S(m,n) = sum_{k=1}^{n-1} {mk/n} cot(k pi/n)   (Vasyunin / BDBLS p.141)"""
    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):
    """A(b,r) for coprime b,r: <f_b,f_r> = [C(b+r) + ((b-r)/2)log(r/b) - (pi/2)(S(b,r)+S(r,b))]/(br)"""
    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)

def mu_sieve(n):
    mu=np.ones(n+1,dtype=np.int64); mu[0]=0
    primes=[]; isc=np.zeros(n+1,dtype=bool)
    for i in range(2,n+1):
        if not isc[i]: primes.append(i); mu[i]=-1
        for p in primes:
            if i*p>n: break
            isc[i*p]=True
            if i%p==0: mu[i*p]=0; break
            else: mu[i*p]=-mu[i]
    return mu

NMAX = 2000
mu = mu_sieve(NMAX)
sq = [n for n in range(1,NMAX+1) if mu[n]!=0]
print(f"squarefree n <= {NMAX}: {len(sq)}", flush=True)

t0=time.time()
idx={n:i for i,n in enumerate(sq)}
K=len(sq)
Gs=np.zeros((K,K))
for i,m in enumerate(sq):
    for j in range(i,K):
        n=sq[j]
        g=gcd(m,n)
        Gs[i,j]=Gs[j,i]=Fexpr(n//g, m//g)/g
    if i%100==0: print(f"  row {i}/{K}  [{time.time()-t0:.0f}s]", flush=True)
print(f"Gram (squarefree) done [{time.time()-t0:.0f}s]", flush=True)
np.save('nb_gram_sq.npy', Gs); json.dump(sq, open('nb_sq.json','w'))

c = np.array([(GAMMA-1-log(n))/n for n in sq])
muv = np.array([mu[n] for n in sq], dtype=float)
logn = np.array([log(n) for n in sq])

res=[]
for N in [10,20,50,100,200,400,700,1000,1500,2000]:
    mask = np.array([n<=N for n in sq])
    a = np.where(mask, muv*(1-logn/log(N)), 0.0)
    E = 1 - 2*float(a@c) + float(a@Gs@a)
    res.append((N, E, E*log(N)))
    print(f"N={N:5d}:  E(V_N)={E:.6f}   E*logN={E*log(N):.5f}   target 2B^={C_BD:.5f}", flush=True)
json.dump(res, open('nb_VN.json','w'))

print("\noptimal d_N^2 (all n <= N):", flush=True)
resd=[]
for N in [5,10,20,40,80,160,320]:
    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)
    wmin=w.min()
    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)
    resd.append((N,d2,d2*log(N),float(wmin)))
    print(f"N={N:4d}: d_N^2={d2:.6f}  d^2*logN={d2*log(N):.5f}  min_eig={wmin:.2e}", flush=True)
json.dump(resd, open('nb_dN.json','w'))
print("DONE", flush=True)