back to table · edit · history · where entries came from · files · download
5132 bytes, as of the version from 2026-10-27 07:39 (current). Recorded here, not run.
"""Numerical checks for T541, separate full-trace formulas and controls.
Run: sage -python checks.py [--stored]
Download generate.py, accelerate.py, and check.py alongside this file.
The high-precision Voros check shares only the generic summation helper with
the Watkins computation. An optional --library check instead uses mpmath's
own hypergeometric implementation at 25 decimal digits on every N>=3.
Neither high-precision agreement nor that library check proves a tail bound.
"""
import os,sys,json,time
import numberdb.sage as numberdb
import mpmath as mp
from sage.rings.real_mpfi import RealIntervalField
from generate import compute,boundary,FullSpectralZeta
from accelerate import hyper_tail
from check import agrees_with,measure,exactness,names_its_rings
# This formula is transcribed independently from Voros 1983 (C.22)-(C.28).
def voros(N,n,working=190,library=False):
with mp.workdps(working):
if N==1:
# Independently take log derivatives of the Airy determinants.
return mp.nstr(-mp.diff(lambda E:mp.log(mp.airyai(-E)*(-mp.airyai(-E,1))),0,2),working)
if N==2:
# Spectral odd-integer Hurwitz-zeta representation.
return mp.nstr(mp.zeta(n,mp.mpf(1)/2)/mp.mpf(2)**n,working)
v=mp.mpf(1)/(N+2);h=mp.mpf(1)/2
parity=mp.sin(mp.pi*v)/(2*mp.sqrt(mp.pi))*(2*v)**(2-4*v)
parity*=mp.gamma(v)*mp.gamma(2*v)*mp.gamma(3*v)/mp.gamma(2*v+h)
if n==1:
return mp.nstr(mp.tan(2*mp.pi*v)/mp.tan(mp.pi*v)*parity,working)
A=[4*v,2*v,3*v,5*v,h]
B=[1+2*v,1+v,1-v,h+4*v]
f=mp.hyper(A,B,1) if library else hyper_tail(A,B,240,160)
pref=mp.sin(mp.pi*v)/(4*v*v*mp.sqrt(mp.pi))*(2*v)**(4-8*v)
pref*=mp.gamma(3*v)*mp.gamma(4*v)*mp.gamma(5*v)/mp.gamma(h+4*v)
return mp.nstr(parity**2+pref*f,working)
def close(a,b,tolerance):
return abs(a-b)<tolerance*max(mp.mpf(1),abs(b))
def run_checks(stored_mode=False,library=False):
if os.environ.get('NUMBERDB_KEY_FROM_STDIN')=='1':
os.environ['NUMBERDB_API_KEY']=sys.stdin.read().strip()
mp.mp.dps=220
gen=FullSpectralZeta()
# Known, nonzero control precedes any interpretation of the method.
with mp.workdps(190):
control=hyper_tail([mp.mpf(1)]*2,[mp.mpf(3)],240,160)
assert abs(control-2)<mp.mpf('1e-150')
assert not close(mp.mpf(2)+mp.mpf('1e-90'),mp.mpf(2),mp.mpf('1e-110'))
print('CONTROL 2F1(1,1;3;1)=2 and deliberate 1e-90 perturbation: PASS',flush=True)
from numberdb._generate import _stored_entries
records=_stored_entries(numberdb.table('T541')) if stored_mode else None
field=RealIntervalField(800)
vals={};independent={};comparison={};reference={};minimum=mp.inf
for p in gen.enumerate():
N,n=p['N'],p['n'];key=f'N={N},n={n}'
actual=mp.mpf(compute(N,n,200))
other=mp.mpf(voros(N,n,25 if library else 190,library))
agreement=-mp.log10(abs((actual-other)/actual)) if actual!=other else mp.inf
assert agreement>(20 if library else 115),(key,'Voros disagreement',str(agreement))
assert actual>0
independent[key]=field(str(other))
vals[key]=gen.value(p,100)
compare_digits=18 if library else 100
comparison[key]=mp.nstr(actual,compare_digits)
reference[key]=mp.nstr(other,compare_digits)
if N>=3:
# Watkins 2.10a,b, derived from the quantum Wronskian, checks both
# parity sectors independently of the generic summation formula.
z1m=mp.mpf(boundary(N,1,1,200));z1p=mp.mpf(boundary(N,1,-1,200))
sig=mp.mpf(2)/(N+2);lam=mp.mpf(1)/2
coeff=lambda a:mp.sin(mp.pi*sig*(lam+a))/mp.sin(mp.pi*sig*lam)
if n==1:
residual=coeff(1)*z1m+coeff(-1)*z1p
else:
zm=mp.mpf(boundary(N,2,1,200));zp=mp.mpf(boundary(N,2,-1,200))
residual=coeff(2)*zm+coeff(-2)*zp+(coeff(1)**2-coeff(2))*(z1p-z1m)**2
assert abs(residual)<mp.mpf('1e-120'),(key,'Wronskian',str(residual))
if stored_mode and not library:
text=records[f'{N},{n}']
from decimal import Decimal
center=Decimal(text); ulp=Decimal(10)**center.as_tuple().exponent
mid=field(str(center));rad=field(str(ulp))
interval=field(mid.lower()-rad.upper(),mid.upper()+rad.upper())
assert interval.lower()<independent[key].lower() and independent[key].upper()<interval.upper(),key
minimum=min(minimum,agreement)
print(json.dumps({'key':key,'independent_agreement_digits':str(agreement),'stored':stored_mode,'library':library, 'wronskian_residual':str(residual) if N>=3 else None}),flush=True)
assert len(vals)==22
assert not exactness(vals)
assert not names_its_rings(__file__.replace('checks.py','generate.py'))
assert not agrees_with(comparison, lambda key: reference[key])
print('MINIMUM_INDEPENDENT_AGREEMENT',mp.nstr(minimum,12),'MEASURE',measure(vals),flush=True)
return vals
if __name__=='__main__':
run_checks('--stored' in sys.argv,'--library' in sys.argv)