back to table · edit · history · where entries came from · files · download
6607 bytes, as of the version from 2026-10-25 13:21 (current). Recorded here, not run.
"""Independent checks for Bessel cross-product zeros -- numberdb.org/T487.
Run with SageMath:
sage -python checks.py values.json # flat records, no network
sage -python checks.py --stored # authenticated draft read-back
The generator is neither imported nor called. All written decimal intervals
are checked using Arb at 896 bits. Indexing is checked by a separately coded
Sturm count using a rational mesh. Every root is also checked by integrating
the transformed radial differential equation in Pruefer phase form with
SciPy DOP853, independent of Bessel evaluations. That ODE test is numerical,
not a replacement for the interval certificate. Its control is the known
free equation at nu=1/2, including a perturbed value that must fail.
"""
import numberdb.sage as numberdb
from sage.rings.rational_field import QQ
from sage.rings.complex_arb import ComplexBallField
from sage.rings.real_arb import RealBallField
from scipy.integrate import solve_ivp
import json
import math
import os
import sys
from collections import defaultdict
from pathlib import Path
BITS = 896
C = ComplexBallField(BITS)
B = RealBallField(BITS)
def rational_decimal(text):
mantissa, _, exponent = str(text).lower().partition('e')
a, dot, b = mantissa.partition('.')
if not dot:
raise ValueError('this table has non-rational roots, so decimals are required')
power = int(exponent or '0')-len(b)
return QQ(int(a+b))*QQ(10)**power, QQ(10)**power
def sign_of(v):
assert v.real().is_finite() and v.imag().is_finite()
assert v.imag().contains_zero()
if v.real() > 0:
return 1
if v.real() < 0:
return -1
raise AssertionError('the sign is not certified')
def determinant(order, ratio, x):
z, w, v = C(x), C(ratio*x), C(order)
return z.bessel_J(v)*w.bessel_Y(v)-z.bessel_Y(v)*w.bessel_J(v)
def count_nodes(order, ratio, endpoint):
# Independent mesh: use h <= 1/(endpoint+1), and verify the exact bound.
intervals = int(((ratio-1)*(endpoint+1)).ceil())+2
h = (ratio-1)/intervals
assert h*h*(endpoint*endpoint+QQ(1)/4) < 9
last, changes = 1, 0
for j in range(1,intervals+1):
s = sign_of(determinant(order,1+j*h,endpoint))
changes += int(s != last)
last = s
return changes
def phase_at_boundary(order, ratio, root):
# u'' + (root^2-(order^2-1/4)/r^2)*u = 0.
# u=R*sin(theta), u'=root*R*cos(theta), u(1)=0, u'(1)>0.
c = order*order-0.25
def ode(r,y):
return [root-c*math.sin(y[0])**2/(root*r*r)]
answer=solve_ivp(ode,(1.0,ratio),[0.0],method='DOP853',
rtol=2e-13,atol=2e-14,max_step=(ratio-1)/80)
assert answer.success and math.isfinite(answer.y[0,-1])
return float(answer.y[0,-1])
def read_stored():
doc = numberdb.table('T487')
records=[]
for nu,rows in doc['Numbers'].items():
for ratio,entries in rows.items():
for k,entry in entries.items():
value=entry['number'] if isinstance(entry,dict) else entry
records.append({'params':{'nu':str(nu),'lambda':str(ratio),'k':str(k)},'number':value})
return records
def run(records):
# Controls precede measurement. Nonzero pi is required, not just zero residual.
control=phase_at_boundary(0.5,2.0,math.pi)
assert abs(control-math.pi) < 1e-12
perturbed=phase_at_boundary(0.5,2.0,1.01*math.pi)
assert abs(perturbed-math.pi) > 0.03
for x,expected in [(QQ(3),0),(QQ(4),1),(QQ(7),2)]:
assert count_nodes(QQ(1)/2,QQ(2),x)==expected
# Analytic half-integer identity is an external control for the Bessel determinant.
for ratio in map(QQ,['5/4','3/2','2','5/2','3','4','5','10']):
x=QQ(7)/3
actual=determinant(QQ(1)/2,ratio,x).real()
known=2*(B(ratio-1)*B(x)).sin()/(B.pi()*B(ratio).sqrt()*B(x))
assert actual.is_finite() and known.is_finite()
assert B(actual.rad()) < B(10)**(-140) and B(known.rad()) < B(10)**(-140)
assert actual.overlaps(known)
print(json.dumps({'controls':'passed','known_phase':control,'perturbed_phase':perturbed}),flush=True)
groups=defaultdict(list)
max_residual=0.0
residual_at=None
for record in records:
params=record['params']
nu,ratio,k=QQ(params['nu']),QQ(params['lambda']),int(params['k'])
x,ulp=rational_decimal(record['number'])
a,b=x-ulp,x+ulp
assert a>0
s1,s2=sign_of(determinant(nu,ratio,a)),sign_of(determinant(nu,ratio,b))
assert s1*s2==-1
phase=phase_at_boundary(float(nu),float(ratio),float(x))
residual=abs(phase-k*math.pi)
assert residual<1e-8, (params,residual)
if residual>max_residual:
max_residual,residual_at=residual,params
c=nu*nu-QQ(1)/4
base=(B(k)*B.pi()/B(ratio-1))**2
qmin=min(c,c/(ratio*ratio))
qmax=max(c,c/(ratio*ratio))
squared=B(x,ulp)**2
assert squared>base+B(qmin) and squared<base+B(qmax)
groups[(str(nu),str(ratio))].append((k,a,b,x,ulp))
assert len(records)==480 and len(groups)==48
for (nu,ratio),rows in groups.items():
rows.sort()
assert [r[0] for r in rows]==list(range(1,11))
assert all(a[2]<b[1] for a,b in zip(rows,rows[1:]))
assert count_nodes(QQ(nu),QQ(ratio),rows[-1][2])==10
# T390 is an independent corpus value of the annulus eigenvalue.
old=numberdb.table('T390')['Numbers']['annulus-1-2']['1']
if isinstance(old,dict):
old=old['number']
target,target_ulp=rational_decimal(old)
first=next(r for r in groups[('0','2')] if r[0]==1)
square=B(first[3],first[4])**2
assert square.is_finite() and B(square.rad())<B(10)**(-95)
assert square.overlaps(B(target,target_ulp))
result={'entries':len(records),'certified_decimal_brackets':len(records),
'minmax_bounds':len(records),'certified_index_groups':len(groups),'independent_ode_checks':len(records),
'max_phase_residual':max_residual,'max_phase_residual_at':residual_at,
'ode_phase_tolerance':1e-8,'T390_annulus_square':'agrees',
'bits':BITS,'controls':'passed'}
print(json.dumps(result),flush=True)
return result
if __name__=='__main__':
if os.environ.get('NUMBERDB_KEY_FROM_STDIN')=='1':
numberdb.configure(api_key=sys.stdin.read().strip())
if '--stored' in sys.argv:
records=read_stored()
else:
path=sys.argv[1] if len(sys.argv)>1 else str(Path(__file__).with_name('values.json'))
records=json.loads(Path(path).read_text())
run(records)