back to table · edit · history · where entries came from · files · download
7615 bytes, as of the version from 2026-11-09 18:56 (current). Recorded here, not run.
"""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()