back to table · edit · history · where entries came from · files · download
8416 bytes, as of the version from 2026-11-09 14:49 (current). Recorded here, not run.
"""Independent coefficient, numerical and identity checks for T690.
Run: sage -python checks.py
For checks on the saved draft: NUMBERDB_CHECK_STORED=1 sage -python checks.py
Set NUMBERDB_API_KEY for private draft reads. Output is newline-delimited JSON
prefixed by CHECK, including the independent values and individual residuals.
The primary generator uses cubic Euler recurrences and PARI lfun. This checker
uses Sage's modular symbols eigenforms and Dokchitser's implementation, whose
symmetric-square coefficients use a different Dirichlet convolution. PARI's
Haberland period pairing checks the Petersson relation for every form.
"""
import numberdb.sage as numberdb
from sage.modular.modform.constructor import Newforms
from sage.rings.real_mpfr import RealField
from sage.rings.rational_field import QQ
from sage.libs.pari import pari
from sage.env import SAGE_VERSION
import json, os, sys, math, pathlib
import generate as g
R=RealField(900)
def emit(kind, **fields):
print('CHECK',json.dumps(dict(kind=kind,**fields)),flush=True)
def relative(x,y):
x,y=R(x),R(y)
assert not (x.is_NaN() or x.is_infinity() or y.is_NaN() or y.is_infinity()) and x!=0 and y!=0
return abs((x-y)/y)
def require_near(x,y,digits,label):
e=relative(x,y)
assert e<R(10)**(-digits),(label,str(e))
return str(e)
def gamma_factor(k,s):
z=R(s)
return 2*(2*R.pi())**(-z)*z.gamma()*R.pi()**(-(z-k+2)/2)*((z-k+2)/2).gamma()
def read_stored():
doc=numberdb.table('T690');block=doc['Numbers'];out={}
if isinstance(block,list):
for row in block:out[tuple(int(row['params'][x]) for x in ('k','i','s'))]=row['number']
else:
for k,byi in block.items():
for i,bys in byi.items():
for s,row in bys.items():out[(int(k),int(i),int(s))]=row['number'] if isinstance(row,dict) else row
return out
def main():
if os.environ.get('NUMBERDB_KEY_FROM_STDIN')=='1':numberdb.configure(api_key=sys.stdin.read().strip())
pari.allocatemem(160*1024*1024,silent=True)
emit('versions',sage=SAGE_VERSION,pari=str(pari('version()')),numberdb=numberdb.__version__)
stored=read_stored() if os.environ.get('NUMBERDB_CHECK_STORED')=='1' else None
if stored is not None:assert set(stored)=={tuple(p[x] for x in ('k','i','s')) for p in g.SymmetricSquareValues().enumerate()}
# Controls must be nonzero known answers, before numerical residuals matter.
assert relative(2,1)==1
assert abs(R(pari.lfun(1,2,precision=700))-R.pi()**2/6)<R(10)**-200
fcontrol=Newforms(1,12,names='a')[0]
assert list(fcontrol.qexp(6))[1:]==[1,-24,252,-1472,4830]
# Sage manual's 300-digit regression value is an external numerical control.
delta22='0.999645711124771397835729622033153189549796658647254961493709341358991830134499267117001769570658192128781135161587571716303826382489492569725002840546129937149159065273765309218543427544527498868033604310899372849565046516553245752253255585377793879866297612679545029546953895098375829822346290125161'
emit('controls',nonzero_relative_residual='1',zeta2='passed',delta_fourier='passed')
roster=[]
for k in g.WEIGHTS:
forms=Newforms(1,k,names='a');assert len(forms)==1
f=forms[0];K=f.base_ring();a2=f[2]
degree=int(K.degree())
assert degree==(2 if k==24 else 1)
nmax=650
a=g.fourier_coefficients(k,nmax)
fq=f.qexp(nmax+1)
assert all(K(A)+K(B)*a2==fq[n] for n,(A,B) in enumerate(a))
assert all(int(pari.mfdim([1,j],1))==(2 if j==24 else (1 if j in g.WEIGHTS else 0)) for j in range(2,27,2))
bound=k//12
# This exact equality identifies our modular form against an independent
# eigenform construction beyond the Sturm bound, not just its name.
embs=sorted(K.embeddings(RealField(700)),key=lambda e:e(a2))
mf=pari.mfinit([1,k],0);pf=pari.mfeigenbasis(mf)[0]
pc=pari.mfcoefs(pf,nmax)
if k==24:
# PARI expresses a_2=-24*y+552, y^2-y-36042=0.
p2=pc[2]
assert p2.minpoly()==pari('x^2-1080*x-20468736')
assert all(pari(A)+pari(B)*p2==pc[n] for n,(A,B) in enumerate(a))
else:
assert all(pari(A)==pc[n] and B==0 for n,(A,B) in enumerate(a))
# Exact identity L(sym^2,s)=zeta(2s-2k+2)/zeta(s-k+1)
# * sum a_n^2/n^s.
b=g.symmetric_coefficients(k,nmax)
squares=[g.mul(v,v) for v in a]
mu=[0]+[int(pari.moebius(n)) for n in range(1,nmax+1)]
conv=[g.ZERO]*(nmax+1)
for d in range(1,nmax+1):
factor=mu[d]*d**(k-1)
for n in range(d,nmax+1,d):conv[n]=g.add(conv[n],g.scale(squares[n//d],factor))
alternative=[g.ZERO]*(nmax+1)
for r in range(1,math.isqrt(nmax)+1):
for n in range(r*r,nmax+1,r*r):alternative[n]=g.add(alternative[n],g.scale(conv[n//(r*r)],r**(2*k-2)))
assert alternative[1:]==list(b[1:])
# Gamma regularity derives the critical set independently of the enumeration.
def regular(s):return s>0 and not (s-k+2<=0 and (s-k+2)%2==0)
assert [s for s in range(-100,101) if regular(s) and regular(2*k-1-s)]==list(g.critical_integers(k))
# Period pairing is independent of both symmetric-square evaluators.
symbols=pari.mfsymbol(mf,pf,precision=700)
norms=pari.mfpetersson(symbols)
if k==24:
roots=list(pari.polrootsreal(pari.mfparams(pf)[3],precision=700))
norm_order=sorted(range(2),key=lambda j:R(552)-24*R(roots[j]))
norm_values=[R(norms[j,j]) for j in norm_order]
else:norm_values=[R(norms)]
for i,emb in enumerate(embs,1):
low,lmeta=g.compute_form(k,i,160)
high,hmeta=g.compute_form(k,i,220)
L=f.symsquare_lseries(embedding=(emb if degree>1 else 0),prec=700)
feq=L.check_functional_equation()
emit('form',k=k,i=i,sturm_bound=bound,exact_coefficients_checked=nmax,
embedded_a2=str(R(emb(a2))),pari160=lmeta,pari220=hmeta,dok_feq=str(feq))
roster.append({'k':k,'i':i,'sturm_bound':bound,
'a2':str(a2) if degree==1 else {'polynomial':'x^2-1080*x-20468736','isolating_interval':[-4017,-4016] if i==1 else [5096,5097]},
'coefficients_0_through_sturm_as_pairs':a[:bound+1],
'pair_meaning':'A+B*a2', 'critical_integers':g.critical_integers(k)})
independent={s:R(L(s).real()) for s in g.critical_integers(k)}
for s in g.critical_integers(k):
x=R(high[s]);y=independent[s]
# 110 digits for the primary precision comparison leaves a measured
# margin over the 100 displayed significant digits.
e_precision=require_near(low[s],x,110,(k,i,s,'precision'))
e_independent=require_near(y,x,150,(k,i,s,'Dokchitser'))
t=2*k-1-s
e_reflection=require_near(gamma_factor(k,s)*x,gamma_factor(k,t)*R(high[t]),170,(k,i,s,'functional equation'))
row={'k':k,'i':i,'s':s,'pari220':high[s],'dok700':str(y),
'relative_precision_difference':e_precision,'relative_independent_difference':e_independent,
'relative_feq_difference':e_reflection}
if stored is not None:
from decimal import Decimal
text=str(stored[(k,i,s)])
value=R(text);ulp=R(10)**Decimal(text).as_tuple().exponent
assert abs(value-y)<=ulp,(k,i,s,'stored ulp')
row['stored_ulp_error']=str(abs(value-y)/ulp)
assert abs(gamma_factor(k,s)*value-gamma_factor(k,t)*R(stored[(k,i,t)]))<=abs(gamma_factor(k,s))*ulp+abs(gamma_factor(k,t))*R(10)**Decimal(str(stored[(k,i,t)])).as_tuple().exponent
emit('entry',**row)
norm=norm_values[i-1]
predicted=R(math.factorial(k-1))*R(high[k])/(R(2)**(2*k-1)*R.pi()**(k+1))
e_norm=require_near(predicted,norm,150,(k,i,'Petersson'))
if k==12:require_near(high[22],delta22,180,'Sage regression delta22')
emit('norm',k=k,i=i,petersson=str(norm),relative_difference=e_norm)
emit('roster',forms=roster)
emit('complete',entries=162,stored_checked=stored is not None)
if __name__=='__main__':main()