back to table · edit · history · where entries came from · files · download
7393 bytes, as of the version from 2026-11-09 18:56 (current). Recorded here, not run.
"""Zeros of period polynomials of newforms -- numberdb.org/T694.
Run with SageMath, numberdb, and the attached roster.json:
sage -pip install numberdb
sage -python generate.py # verify every stored root
sage -python generate.py --publish # fill the existing draft
All 626 roots of all 121 embedded forms in the four slices described in
argument-selection.md are enumerated. Exact Sturm coefficients and isolating
intervals identify forms independently of PARI's list order. The polynomial
is integral_0^(i infinity) f(tau)(tau-X)^(k-2) d tau, without rescaling X.
PARI mfperiodpol uses extra working precision, 160 and 220 decimal digits.
Componentwise agreement supports 100 digits; it is not a certified bound.
checks.py independently reconstructs its coefficients by split Mellin sums.
"""
import sys,json,math
from pathlib import Path
from functools import lru_cache
from fractions import Fraction
import numberdb.sage as numberdb
from numberdb import ComplexInterval, RealInterval
from sage.libs.pari import pari
from sage.rings.rational_field import QQ
from sage.rings.integer_ring import ZZ
from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
from sage.rings.qqbar import AA, _init_qqbar
from sage.rings.complex_mpfr import ComplexField
from sage.rings.real_mpfr import RealField
from sage.rings.complex_interval_field import ComplexIntervalField
from sage.rings.real_arb import RealBallField
from sage.rings.real_mpfi import RealIntervalField
_init_qqbar()
BASE=Path(__file__).resolve().parent
ROSTER=json.loads((BASE/'roster.json').read_text())
BY_KEY={(r['N'],r['k'],r['i']):r for r in ROSTER}
P=PolynomialRing(QQ,'x')
WORKING_DIGITS=(160,220)
def coefficient(c):return P(str(pari.lift(c)).replace('y','x'))
@lru_cache(None)
def space(N,k):return pari.mfinit([N,k],0)
@lru_cache(None)
def reconstruct(N,k,i):
row=BY_KEY[N,k,i];F=P(row['field_poly']);cs=[P([QQ(c) for c in v]) for v in row['coefficients']]
lo,hi=map(QQ,row['root_interval']);roots=[z for z in F.roots(AA,multiplicities=False) if lo<z<hi]
assert len(roots)==1
mf=space(N,k);found=[]
for o,f in enumerate(pari.mfeigenbasis(mf)):
if coefficient(pari.mfparams(f)[3])==F and [coefficient(c) for c in pari.mfcoefs(f,row['sturm_bound'])]==cs:found.append((o,f))
assert len(found)==1,(N,k,i)
o,f=found[0]
A=pari.mfatkininit(mf,N,precision=128);assert A[2]==1
assert all(a==row['fricke']*b for a,b in zip(pari.mfcoefs(pari.mfatkin(A,f),row['sturm_bound']),pari.mfcoefs(f,row['sturm_bound'])))
return mf,f,o,roots[0]
@lru_cache(None)
def orbit_polynomials(N,k,o,working):
bits=numberdb.bits(working,losing=32);pari.set_real_precision(working+12)
mf=space(N,k);f=pari.mfeigenbasis(mf)[o];d=int(pari.poldegree(pari.mfparams(f)[3]))
raw=pari.mfperiodpol(mf,f,precision=bits)
return [raw] if d==1 else list(raw)
def embedding_position(f,z,d,bits):
if d==1:return 0
C=ComplexField(bits);images=[C(v) for v in pari.mfembed(f,pari.Mod('y',pari.mfparams(f)[3]),precision=bits)]
hits=[j for j,v in enumerate(images) if abs(v-C(z))<C(2).real()**(-bits+60)]
assert len(hits)==1
return hits[0]
def ordered_roots(coeff,row,working):
bits=numberdb.bits(working,losing=32);C=ComplexField(bits);R=RealField(bits)
poly=pari.Polrev([pari(c) for c in coeff]);roots=[C(z) for z in pari.polroots(poly,precision=bits)]
assert len(roots)==row['k']-2
# These replacements follow identities; proximity checks only locate the
# already proved roots among the numerical output, never prove exactness.
forced=[]
if row['N']==1:forced += [(C(1),'real+'),(C(-1),'real-')]
if row['N']==1 and (row['k']-2)%6==0:
for re in [-1,1]:
for im in [-1,1]:forced.append((C(R(re)/2,R(im)*R(3).sqrt()/2),'elliptic:%s:%s'%(re,im)))
if row['epsilon']==-1:
y=1/R(row['N']).sqrt();forced += [(C(0,y),'imag+'),(C(0,-y),'imag-')]
labels={};tol=R(10)**(-working+35)
for target,kind in forced:
hits=[j for j,z in enumerate(roots) if abs(z-target)<tol];assert len(hits)==1,(row,kind)
j=hits[0];assert j not in labels;labels[j]=kind;roots[j]=target
for j,z in enumerate(roots):
assert not any(t.is_NaN() or t.is_infinity() for t in [z.real(),z.imag()])
if j not in labels:
assert abs(z.real())>tol and abs(z.imag())>tol,('unproved axis root',row['N'],row['k'],row['i'],str(z))
def angle(z):
a=z.argument();return a if a>=0 else a+2*R.pi()
out=sorted([(z,labels.get(j)) for j,z in enumerate(roots)],key=lambda v:(angle(v[0]),abs(v[0])))
assert min(abs(a[0]-b[0]) for j,a in enumerate(out) for b in out[:j])>R('1e-50')
return out
@lru_cache(None)
def calculation(N,k,i,working):
row=BY_KEY[N,k,i];bits=numberdb.bits(working,losing=32);C=ComplexField(bits)
pari.set_real_precision(working+12)
mf,f,o,z=reconstruct(N,k,i);j=embedding_position(f,z,row['field_degree'],bits)
pp=orbit_polynomials(N,k,o,working)[j]
coeff=[C(pari.polcoef(pp,n)) for n in range(k-1)]
roots=ordered_roots(coeff,row,working)
return {'coefficients':coeff,'roots':roots}
def interval_part(ball):
interval=RealIntervalField(ball.parent().precision())(ball)
def exact(t):
q=t.exact_rational();return Fraction(int(q.numerator()),int(q.denominator()))
return RealInterval(exact(interval.lower()),exact(interval.upper()))
class PeriodPolynomialZeros(numberdb.Generator):
table='T694'
parameters=('N','k','i','j')
type='C'
digits=100
rigour='heuristic (agreement-checked)'
files=('generate.py','roster.json','argument-selection.md','checks.py','prepare.py','companion-roster.json','dimensions.json','numerical-evidence.json','skill-evidence.json','check.py','dry_run.py','independent-roots.json','polynomials.json')
def enumerate(self):
for r in ROSTER:
for j in range(1,r['k']-1):yield {'N':r['N'],'k':r['k'],'i':r['i'],'j':j}
def value(self,params,digits):
assert digits<=100
N,k,i,j=(int(params[x]) for x in self.parameters);row=BY_KEY[N,k,i]
vals={w:calculation(N,k,i,w)['roots'][j-1] for w in WORKING_DIGITS}
kind=vals[160][1];assert vals[220][1]==kind
if kind in ['real+','real-']:
return {'number':ZZ(1 if kind=='real+' else -1),'comment':'This exact root follows from the level-one period relations.'}
if kind and kind.startswith('elliptic:'):
unused,re,im=kind.split(':');R=RealBallField(numberdb.bits(130))
v=ComplexInterval(Fraction(int(re),2),interval_part(int(im)*R(3).sqrt()/2))
return {'number':v,'comment':r'The order-three level-one period relation forces this root $(%s%s\mathrm{i}\sqrt3)/2$.'%(re,'+' if im=='1' else '-')}
if kind in ['imag+','imag-']:
sign=1 if kind=='imag+' else -1
if math.isqrt(N)**2==N:v=ComplexInterval(0,Fraction(sign,math.isqrt(N)))
else:
R=RealBallField(numberdb.bits(130));v=ComplexInterval(0,interval_part(sign/R(N).sqrt()))
return {'number':v,'comment':r'The functional-equation sign is $-1$, which forces the root $'+('i' if sign==1 else '-i')+r'/\sqrt{'+str(N)+r'}$.'}
real=numberdb.agreeing(lambda w:vals[w][0].real(),at=WORKING_DIGITS)
imag=numberdb.agreeing(lambda w:vals[w][0].imag(),at=WORKING_DIGITS)
return ComplexIntervalField(numberdb.bits(220))(real,imag)
if __name__=='__main__':
if '--publish' in sys.argv:print(PeriodPolynomialZeros().publish(message='Compute all period-polynomial zeros with exact form identities and independent Mellin checks'))
else:
report=PeriodPolynomialZeros().verify(sample=None);print(report);sys.exit(0 if report.ok else 1)