import numpy as np, json, sys
GAPFREE=len(sys.argv)>1
M=1<<16; th=2*np.pi*(np.arange(M)+0.5)/M
f=lambda p: np.log(np.abs(2*np.sin(p/2)))
def Tm(vals,n):
    c=np.fft.fft(vals)/M  # c[k] ~ coefficient of e^{-ik th}... use symmetric indexing
    idx=np.arange(n); K=idx[:,None]-idx[None,:]; return c[K % M]
def c_of_s(N,s):
    n=N-2; d=2*np.pi*s/N
    w=np.abs(np.exp(1j*th)-np.exp(1j*d/2))**2*np.abs(np.exp(1j*th)-np.exp(-1j*d/2))**2
    g=f(-d/2-th)-f(d/2-th)
    if GAPFREE: w=w*((th>d/2)&(th<2*np.pi-d/2))
    T0=Tm(w,n); T1=Tm(w*g,n); T2=Tm(w*g*g,n); Ti=np.linalg.inv(T0)
    A=Ti@T1; E=np.trace(Ti@T2)+np.trace(A)**2-np.trace(A@A)
    return float(np.real(E))/s**2, float(np.real(np.trace(A)))
out={}
for N in [4,12,25,46]:
    ss=[1e-3,0.05,0.15,0.25,0.35,0.45,0.55,0.75,0.95]; out[N]=[(s,)+c_of_s(N,s) for s in ss]
    print(N,[(s,round(c,4)) for s,c,m in out[N]],'pi^2/15(1-4/N^2)=',round(np.pi**2/15*(1-4/N**2),4))
json.dump(out,open('c_exact_s'+('_gapfree' if GAPFREE else '')+'.json','w'),indent=1)
# simulation comparison: E[(x-y)^2]/s^2 in bins
for N in [12,25,46]:
    a=np.load(f'cue_N{N}.npy'); s,l1,l2=a.T; Y=(l1-l2)**2/s**2
    print('sim',N,[(lo,round(Y[(s>=lo)&(s<lo+.1)].mean(),3),round(Y[(s>=lo)&(s<lo+.1)].std()/np.sqrt(((s>=lo)&(s<lo+.1)).sum()),3)) for lo in [0,.1,.2,.3,.4,.5,.7,.9]])