"""Independent split-Mellin and root checks for numberdb.org/T694.
Run: sage -python checks.py with generate.py, roster.json, check.py present.
The independent analytic path does not call mfperiodpol or lfun. It obtains
exact Fourier coefficients of the identified form, then sums integer
incomplete-gamma expansions at 150 digits/450 terms and 190 digits/600 terms.
Newton refinement uses these separately assembled coefficients. Its starting
points are the main roots; distinctness and the degree account for every root.
These finite sums are agreement tests, not certified tail estimates.
"""
import json,math,sys
from pathlib import Path
from functools import lru_cache
import numberdb.sage as numberdb
from sage.libs.pari import pari
from sage.rings.rational_field import QQ
from sage.rings.real_mpfr import RealField
from sage.rings.complex_mpfr import ComplexField
from sage.env import SAGE_VERSION
from generate import ROSTER,reconstruct,calculation,PeriodPolynomialZeros
from check import exactness,measure,names_its_rings

@lru_cache(None)
def exact_fourier(N,k,o,count):
 mf=pari.mfinit([N,k],0);f=pari.mfeigenbasis(mf)[o]
 return list(pari.mfcoefs(f,count))

def eta_control(row,factors):
 a=[0]*101;a[0]=1
 for scale,power in factors:
  for n in range(1,100//scale+1):
   for unused in range(power):
    for j in range(100,scale*n-1,-1):a[j]-=a[j-scale*n]
 shift=sum(s*p for s,p in factors)//24;assert shift==1
 known=[0]+a[:100]
 mf,f,o,z=reconstruct(row['N'],row['k'],row['i'])
 assert list(pari.mfcoefs(f,100))==known,(row['N'],row['k'],[(n,str(v),known[n]) for n,v in enumerate(pari.mfcoefs(f,100)) if v!=known[n]][:2])
 return 100

@lru_cache(None)
def weights(N,k,working,count):
 R=RealField(numberdb.bits(working,losing=32));T=1/R(N).sqrt();pi=R.pi()
 out={s:[] for s in range(1,k)}
 for n in range(1,count+1):
  c=2*pi*n;x=c*T;ex=(-x).exp();term=R(1);total=term
  for s in range(1,k):
   if s>1:term*=x/R(s-1);total+=term
   out[s].append(R(math.factorial(s-1))*ex*total/c**s)
 # Known nonzero exponential integral controls the first weight.
 assert abs(out[1][0]-(-2*pi*T).exp()/(2*pi))<R(10)**(-working)
 assert out[1][0]>0
 return out

def independent(row,working=190,count=600):
 N,k,i=row['N'],row['k'],row['i'];mf,f,o,z=reconstruct(N,k,i)
 R=RealField(numberdb.bits(working,losing=32));C=ComplexField(R.precision());root=R(z)
 # Evaluate exact coefficient polynomials with a separate Horner loop.
 a=[]
 for v in exact_fourier(N,k,o,count)[1:]:
  vv=pari.lift(v);degree=0 if vv==0 else max(0,int(pari.poldegree(vv)));cs=[QQ(pari.polcoef(vv,j)) for j in range(degree+1)]
  result=R(0)
  for c in reversed(cs):result=result*root+R(c)
  a.append(result)
 W=weights(N,k,working,count)
 I={s:sum((aa*w for aa,w in zip(a,W[s])),R(0)) for s in W}
 M={s:I[s]+row['epsilon']*R(N)**(k//2-s)*I[k-s] for s in I}
 cs=[C(0)]*(k-1)
 for j in range(k-1):cs[k-2-j]=math.comb(k-2,j)*C(0,1)**(1-j)*M[j+1]
 return cs,M

def horner(cs,z):
 p=cs[-1];d=z.parent()(0)
 for c in reversed(cs[:-1]):d=d*z+p;p=p*z+c
 return p,d

def refine(cs,z):
 for step in range(12):
  p,d=horner(cs,z);assert abs(d)>0
  dz=p/d;z-=dz
  if abs(dz)<z.parent()(10).real()**(-170):break
 p,d=horner(cs,z);assert abs(p/d)<z.parent()(10).real()**(-165)
 return z

def diff(a,b):return abs(a-b)/max(1,abs(a),abs(b))

def main():
 controls={}
 for key,factors in [((1,12,1),[(1,24)]),((5,4,1),[(1,4),(5,4)]),((8,4,1),[(2,4),(4,4)]),((3,6,1),[(1,6),(3,6)])]:
  row=next(r for r in ROSTER if (r['N'],r['k'],r['i'])==key)
  controls[str(key)]=eta_control(row,factors)
 C=ComplexField(720);assert abs(refine([C(-2),C(0),C(1)],C(1))-C(2).sqrt())<C(10).real()**(-165)
 print('CONTROLS',controls,flush=True)
 maxima={name:0 for name in ['precision_roots','independent_roots','coefficient_agreement','mellin_truncation','circle','reflection','functional_equation','vieta','manin']}
 details=[];entries={};root_evidence=[];all_polynomials=[];minsep=1
 for row in ROSTER:
  N,k,i=row['N'],row['k'],row['i'];lo=calculation(N,k,i,160);hi=calculation(N,k,i,220)
  small,Msmall=independent(row,150,450);other,M=independent(row,190,600);C=other[0].parent();R=RealField(C.precision());w=k-2
  for a,b in zip(hi['coefficients'],other):maxima['coefficient_agreement']=max(maxima['coefficient_agreement'],diff(a,b))
  for a,b in zip(small,other):maxima['mellin_truncation']=max(maxima['mellin_truncation'],diff(a,b))
  roots=[z for z,kind in hi['roots']];cs=hi['coefficients']
  if N==1:
   relation=list(cs)
   for n,c in enumerate(cs):
    for a in range(n+1):relation[w-n+a]+=c*math.comb(n,a)*(-1)**(n-a)
    for a in range(w-n+1):relation[a]+=c*(-1)**n*math.comb(w-n,a)*(-1)**(w-n-a)
   maxima['manin']=max(maxima['manin'],max(abs(c) for c in relation)/max(1,max(abs(c) for c in cs)))
  sep=min(abs(a-b) for j,a in enumerate(roots) for b in roots[:j]);minsep=min(minsep,sep)
  prod=[C(1)]
  for z in roots:
   new=[C(0)]*(len(prod)+1)
   for t,a in enumerate(prod):new[t]-=z*a;new[t+1]+=a
   prod=new
  for n,c in enumerate(cs):
   maxima['vieta']=max(maxima['vieta'],diff(prod[n],c/cs[-1]))
   fe=-row['epsilon']*C(0,1)**k*(-1)**(w-n)*R(N)**(n-w//2)*cs[w-n]
   maxima['functional_equation']=max(maxima['functional_equation'],diff(c,fe))
  rowerr=0
  for j,((a,kind),(b,kind2)) in enumerate(zip(lo['roots'],hi['roots']),1):
   assert kind==kind2
   z=refine(other,C(b));err=abs(z-b);rowerr=max(rowerr,err)
   maxima['precision_roots']=max(maxima['precision_roots'],abs(a-b))
   maxima['independent_roots']=max(maxima['independent_roots'],err)
   maxima['circle']=max(maxima['circle'],abs(N*abs(b)**2-1))
   maxima['reflection']=max(maxima['reflection'],min(abs(t+b.conjugate()) for t in roots))
   params={'N':N,'k':k,'i':i,'j':j};value=PeriodPolynomialZeros().value(params,100);entries[str(params)]=value
   root_evidence.append({'params':params,'real':str(z.real()),'imag':str(z.imag()),'forced':kind})
  details.append({'N':N,'k':k,'i':i,'max_root_difference':str(rowerr),'minimum_separation':str(sep)})
  all_polynomials.append({'N':N,'k':k,'i':i,'coefficients':[[str(z.real()),str(z.imag())] for z in cs],'mellin_moments':{s:str(v) for s,v in M.items()}})
  print('CHECKED',N,k,i,'difference',str(rowerr)[:24],flush=True)
 # Published controls for the period-polynomial convention.
 delta=calculation(1,12,1,220)['coefficients'];R=RealField(740)
 even={0:-QQ(36)/691,2:1,4:-3,6:3,8:-1,10:QQ(36)/691}
 odd={1:4,3:-25,5:42,7:-25,9:4}
 for n,c in even.items():assert abs(delta[n]/delta[10]-c/even[10])<R('1e-170')
 for n,c in odd.items():assert abs(delta[n]/delta[9]-QQ(c)/4)<R('1e-170')
 r8=[z for z,kind in calculation(8,4,1,220)['roots']]
 assert min(abs(z-ComplexField(700)('0.17037672','0.30979311')) for z in r8)<R('2e-8')
 assert len(entries)==626
 assert not exactness(entries),exactness(entries)[:3]
 assert not names_its_rings(str(Path(__file__).parent/'generate.py'))
 for name,v in maxima.items():assert v<R('1e-120'),(name,str(v))
 evidence={'table':'T694','entries':626,'forms':121,'working_digits':[160,220],'independent_settings':[[150,450],[190,600]],'max_errors':{k:str(v) for k,v in maxima.items()},'minimum_root_separation':str(minsep),'measurement':measure(entries),'controls':controls,'per_form':details,'sage_version':SAGE_VERSION,'pari_version':str(pari('version()')),'skill_manifest':'agents/runs/20261109T165350Z-build-skills/manifest.json','skill_sha256':'886390324dd1971e381a43cf3b46b822ab0bc67d62f31d8481e29ec79556746b'}
 print('EVIDENCE '+json.dumps(evidence),flush=True)
 print('ROOTS '+json.dumps(root_evidence),flush=True)
 print('POLYNOMIALS '+json.dumps(all_polynomials),flush=True)
if __name__=='__main__':main()