back to table · edit · history · where entries came from · files · download
5527 bytes, as of the version from 2026-10-26 13:28 (current). Recorded here, not run.
"""Independent numerical checks for Scorer Gi zeros -- numberdb.org/T516.
Run: sage -python checks.py [--stored]
Requires generate.py in this directory. --stored authenticates a draft using
NUMBERDB_API_KEY, or stdin when NUMBERDB_KEY_FROM_STDIN=1. Nothing is written.
The independent route uses certified Arb quadrature in DLMF 9.12.4, replacing
integral_x^infinity Ai by 1/3-integral_0^x Ai (DLMF 9.10.11). It shares neither
the coefficient recurrence nor the root finder with the enclosure calculation.
mpmath at 150 and 200 digits provides two further numerical root checks.
"""
import os
import sys
import json
import platform
import time
import numberdb.sage as numberdb
from sage.rings.complex_arb import ComplexBallField
from sage.rings.rational_field import QQ
from sage.env import SAGE_VERSION
import mpmath as mp
import generate as gen
C = ComplexBallField(700)
RF = gen.RF
TOL = RF(10)**-155
def integral(f,a,b):
z=C.integral(f,a,b,abs_tol=TOL,rel_tol=TOL)
assert z.real().is_finite() and z.imag().is_finite() and z.imag().contains_zero()
r=z.real()
assert RF(r.rad()) < RF(10)**-145
return r
def airy_gi(x):
a=integral(lambda t,_:t.airy_ai(),0,x)
b=integral(lambda t,_:t.airy_bi(),0,x)
z=C(x)
value=(z.airy_bi()*(RF(QQ(1)/3)-a)+z.airy_ai()*b).real()
assert value.is_finite() and RF(value.rad()) < RF(10)**-140
return value
def controls():
poly=integral(lambda t,_:t*t,0,1)
assert (QQ(1)/3 in poly) and poly > 0
# Integrating Ai' must give Ai(-1)-Ai(0), with nonzero answer.
ai=integral(lambda t,_:t.airy_ai_prime(),0,-1)
exact=(C(-1).airy_ai()-C(0).airy_ai()).real()
assert ai.overlaps(exact) and ai > 0
initial=airy_gi(QQ(0))
assert initial.overlaps(gen.gi(0)) and initial > 0
# Check the separately evaluated first and second initial derivatives.
assert gen.gi(0,1).overlaps((C(0).airy_bi_prime()/3).real())
assert (gen.gi(0,2)+1/RF.pi()).contains_zero()
print('CONTROLS: integral t^2=1/3; integral Ai prime; Gi(0), Gi prime(0), Gi second(0)',flush=True)
def decimal_interval(s):
# Stored real decimals mean plus or minus one unit of the last place.
assert isinstance(s,str) and '[' not in s and '+/-' not in s
mantissa,_,exponent=s.lower().partition('e')
places=len(mantissa.partition('.')[2])
unit=QQ(10)**(int(exponent or 0)-places)
q=gen.rational_decimal(s)
return q-unit,q+unit
def run(stored=False):
started=time.monotonic()
controls()
brackets,stats=gen.certify()
rows=None
if stored:
d=numberdb.table('T516')
rows=d['Numbers']
assert len(rows)==12
if isinstance(rows, dict):
assert [int(k) for k in rows] == list(range(1,13))
written=[v['number'] if isinstance(v,dict) else v for v in rows.values()]
else:
assert [int(r['params']['n']) for r in rows]==list(range(1,13))
written=[r['number'] for r in rows]
brackets=[decimal_interval(s) for s in written]
results=[]
for n,(lo,hi) in enumerate(brackets,1):
mid=(lo+hi)/2
left,right=gen.gi(lo),gen.gi(hi)
assert gen.sign(left)*gen.sign(right)==-1
assert gen.sign(gen.local_enclosures(lo,hi)[1])
# Independent quadrature at both exact bracket endpoints.
qleft,qright=airy_gi(lo),airy_gi(hi)
assert gen.sign(qleft)*gen.sign(qright)==-1
assert gen.sign(qleft)==gen.sign(left)
assert gen.sign(qright)==gen.sign(right)
assert qleft.overlaps(left) and qright.overlaps(right)
root_results=[]
for dps in (150,200):
with mp.workdps(dps):
m=mp.mpf(str(mid.numerator()))/int(mid.denominator())
r=mp.findroot(mp.scorergi,(m-mp.mpf('0.01'),m+mp.mpf('0.01')))
rs=mp.nstr(r,dps-10)
rq=gen.rational_decimal(rs)
assert lo<rq<hi
# These numerical checks are distinct from certification.
residual=mp.diff(mp.scorergi,m,2)-m*mp.scorergi(m)+1/mp.pi
assert abs(residual)<mp.mpf(10)**(-dps+15)
identity=mp.scorergi(m)+mp.scorerhi(m)-mp.airybi(m)
assert abs(identity)<mp.mpf(10)**(-dps+15)
root_results.append(rs)
interval=RF(mid).add_error(RF((hi-lo)/2).upper())
ode=gen.gi(interval,2)-interval*gen.gi(interval)+1/RF.pi()
assert ode.is_finite() and ode.contains_zero() and RF(ode.rad()) < RF(10)**-70
item={'n':n,'lo':str(lo),'hi':str(hi),'root_150':root_results[0],
'root_200':root_results[1],
'quadrature_left':str(qleft),'quadrature_right':str(qright),
'series_left_radius':str(left.rad()),
'root_relative_bits':int(interval.accuracy()),
'ode_radius':str(ode.rad())}
results.append(item)
print('PASSED',n,root_results[1][:32], 'bits',item['root_relative_bits'],flush=True)
print('EVIDENCE_JSON '+json.dumps({'table':'T516','stored':stored,'root_count':stats,
'versions':{'sage':SAGE_VERSION,'mpmath':mp.__version__,
'numberdb':numberdb.__version__,'python':platform.python_version()},
'seconds':time.monotonic()-started,'results':results}),flush=True)
if __name__=='__main__':
if os.environ.get('NUMBERDB_KEY_FROM_STDIN')=='1':
os.environ['NUMBERDB_API_KEY']=sys.stdin.read().strip()
run('--stored' in sys.argv or os.environ.get('NUMBERDB_CHECK_STORED')=='1')