back to table · edit · history · where entries came from · files · download
8588 bytes, as of the version from 2026-10-27 01:12 (current). Recorded here, not run.
"""Independent checks for numberdb.org/T512, the regular Coulomb radial zeros.
sage -python checks.py
sage -python checks.py --stored # authenticated API read for a draft
The generator supplies brackets. These checks use separately normalized full
Coulomb functions, analytic derivatives, Tricomi U, and mpmath at 150 and 200
decimal digits. Output is JSON evidence after the EVIDENCE_JSON marker.
"""
import os
import sys
import json
import math
import time
import importlib.metadata
import numberdb.sage as numberdb
from numberdb._write import to_text
from sage.rings.rational_field import QQ
from sage.rings.complex_arb import ComplexBallField
from sage.rings.real_arb import RealBallField
import mpmath as mp
from generate import (CoulombFZeros, scan, regular_scaled, certified_root,
mpq, rational_decimal, written_interval, WORK_BITS)
C = ComplexBallField(WORK_BITS)
R = RealBallField(WORK_BITS)
I = C(0, 1)
TOL = QQ(1)/10**130
def real(z):
assert z.real().is_finite() and z.imag().is_finite()
assert z.imag().contains_zero()
return z.real()
def zero(z):
z = C(z)
assert z.real().is_finite() and z.imag().is_finite()
assert z.real().contains_zero() and z.imag().contains_zero(), str(z)
assert abs(z).upper() < TOL, str(z)
def normalization(ell, eta):
e = R(eta)
base = R(1) if eta == 0 else 2*R.pi()*e/((2*R.pi()*e).exp()-1)
for j in range(1, ell+1):
base *= e*e+j*j
return R(2)**ell*base.sqrt()/math.factorial(2*ell+1)
def full_f(ell, eta, rho, derivatives=False):
r, e = C(rho), C(eta)
p = ell+1
a, b, z = p-I*e, C(2*p), 2*I*r
pref = C(normalization(ell, eta))*r**p*(-I*r).exp()
m0 = z.hypergeometric([a], [b])
f = real(pref*m0)
if not derivatives:
return f
m1 = 2*I*a/b*z.hypergeometric([a+1], [b+1])
m2 = -4*a*(a+1)/(b*(b+1))*z.hypergeometric([a+2], [b+2])
u = p/r-I
fp = real(pref*(u*m0+m1))
fpp = real(pref*((u*u-p/(r*r))*m0+2*u*m1+m2))
return f, fp, fpp
def full_g(ell, eta, rho):
r, e = C(rho), C(eta)
p, z = ell+1, -2*I*r
a, b = p+I*e, C(2*p)
gamma = a.gamma()
pref = (-I)**ell*(C.pi()*e/2).exp()*gamma/abs(gamma)
pref *= (I*r).exp()*z**p
u = z.hypergeometric_U(a, b)
h = pref*u
hp = pref*((p/r+I)*u+2*I*a*z.hypergeometric_U(a+1,b+1))
for v in (h, hp):
assert v.real().is_finite() and v.imag().is_finite()
assert v.real().rad() < TOL and v.imag().rad() < TOL
return h.real(), hp.real(), h.imag()
def identities(ell, eta, rho):
r, e = R(rho), R(eta)
f, fp, fpp = full_f(ell, eta, rho, True)
g, gp, hf = full_g(ell, eta, rho)
zero(hf-f)
zero(fpp+(1-2*e/r-ell*(ell+1)/(r*r))*f)
zero(g*fp-f*gp-1)
rp = (1+e*e/(ell+1)**2).sqrt()
sp = (ell+1)/r+e/(ell+1)
fplus = full_f(ell+1,eta,rho)
zero(fp-sp*f+rp*fplus)
if ell:
rm = (1+e*e/ell**2).sqrt()
sm = ell/r+e/ell
fminus = full_f(ell-1,eta,rho)
zero(fp-rm*fminus+sm*f)
zero(rm*fminus-(sm+sp)*f+rp*fplus)
return max(f.rad(),fp.rad(),fpp.rad(),g.rad(),gp.rad())
def controls():
# Nonzero known answers first; then checks that reject deliberate errors.
f, fp, fpp = full_f(0, QQ(0), QQ(1), True)
zero(f-R(1).sin()); zero(fp-R(1).cos()); zero(fpp+R(1).sin())
g,gp,hf = full_g(0, QQ(0), QQ(1))
zero(g-R(1).cos()); zero(gp+R(1).sin()); zero(g*fp-f*gp-1)
assert f > QQ(4)/5 and fp > QQ(1)/2
try:
zero(-g*fp+f*gp-1)
except AssertionError:
pass
else:
raise AssertionError('Reversed-sign Wronskian control was accepted')
assert written_interval('3.14') == (QQ(313)/100,QQ(315)/100)
# Exact inequalities behind the uniform index proof.
assert QQ(3)/8/(1-QQ(3)/8) == QQ(3)/5
assert QQ(4)/3+QQ(1)/9 == QQ(13)/9
assert QQ(13)/54 < 1
assert R.pi()/R(33).sqrt() > QQ(1)/2
# Published mpmath documentation examples (A&S comparison), ell=0..3.
published = [('1.143337392','0.2937960375'),
('0.9935056752','0.3708676452'),
('0.7212774133','0.3952408216'),
('0.4313159311','0.320534053')]
for ell,(fs,fps) in enumerate(published):
f,fp,_ = full_f(ell,QQ(2),QQ(5),True)
for value,s in ((f,fs),(fp,fps)):
lo,hi=written_interval(s)
assert R(lo) < value < R(hi)
# Product and gamma normalizations must agree for every selected pair.
for ell in range(4):
for eta_s in ('0','1/2','-1/2','1','-1','2','-2'):
eta=QQ(eta_s)
gamma_norm=C(2)**ell*(-C.pi()*C(eta)/2).exp()*abs((ell+1+I*C(eta)).gamma())/math.factorial(2*ell+1)
zero(normalization(ell,eta)-gamma_norm)
# Independent Bessel and sine function identities at all selected couplings eta=0.
for ell in range(4):
for rho in (QQ(1)/2,QQ(1),QQ(2),QQ(5)):
expected=(C.pi()*C(rho)/2).sqrt()*C(rho).bessel_J(C(QQ(2*ell+1)/2))
zero(full_f(ell,QQ(0),rho)-expected)
return 'nonzero sine/cosine, unit Wronskian, rejected reversed sign, decimal interval, exact index inequalities, 8 published values/derivatives, 16 Bessel function controls'
def stored_rows():
d=numberdb.table('T512')
rows=[]
for ell,es in d['Numbers'].items():
for eta,ns in es.items():
for n,v in ns.items():
rows.append(({'ell':str(ell),'eta':str(eta),'n':str(n)},v['number'] if isinstance(v,dict) else v))
assert len(rows)==336
return rows
def main(stored=False):
started=time.monotonic()
control=controls()
print('Controls passed:',control,flush=True)
generator=CoulombFZeros()
rows=[]
if stored:
rows=stored_rows()
else:
for params in generator.enumerate():
value=generator.value(params,100)
if isinstance(value,dict):value=value['number']
rows.append((params,to_text(value,100)))
evidence=[]
maxradius=QQ(0)
for idx,(p,text) in enumerate(rows):
ell,eta,n=int(p['ell']),QQ(p['eta']),int(p['n'])
lo,hi=written_interval(text)
a,b=scan(ell,str(eta))[n-1]
assert a < lo < hi < b
fl,fh=full_f(ell,eta,lo),full_f(ell,eta,hi)
assert fl*(-1)**(n-1)>0 and fh*(-1)**n>0
# Two separately executed precision levels; the second starts at the
# written midpoint, rather than taking the first solver's result.
agreements=[]
for dps in (150,200):
with mp.workdps(dps):
mid=mpq((lo+hi)/2)
answer=mp.findroot(lambda x:mp.coulombf(ell,mpq(eta),x),
(mid-mp.mpf('0.01'),mid+mp.mpf('0.01')),
tol=mp.mpf(10)**(-dps+10))
assert mpq(lo)<answer<mpq(hi)
agreements.append(mp.nstr(answer,130))
if eta == 0:
expected=mp.besseljzero(mp.mpf(ell)+mp.mpf(1)/2,n)
assert mpq(lo)<expected<mpq(hi)
if ell==0:assert mpq(lo)<n*mp.pi<mpq(hi)
assert abs(rational_decimal(agreements[0])-rational_decimal(agreements[1]))<QQ(1)/10**125
radius=identities(ell,eta,(lo+hi)/2)
maxradius=max(maxradius,radius)
evidence.append({'params':p,'number':text,'scan_bracket':[str(a),str(b)],
'mpmath_150':agreements[0],'mpmath_200':agreements[1],
'certified_endpoint_signs':[(-1)**(n-1),(-1)**n]})
if (idx+1)%12==0:
print('Checked',idx+1,'roots; ell,eta =',ell,str(eta),'elapsed',round(time.monotonic()-started,1),flush=True)
summary={'table':'T512','mode':'stored' if stored else 'preflight','entries':len(evidence),
'sign_brackets':len(evidence),'indexed_scans':28,'bessel_zero_checks':48,
'mpmath_precision_checks':672,'ode_wronskian_derivative_checks':336,
'order_recurrence_checks':252,'maximum_function_ball_radius':str(maxradius),
'controls':control,'working_bits':WORK_BITS,'scan_bits':160,
'root_enclosure_width':'2e-120','elapsed_seconds':round(time.monotonic()-started,1),
'versions':{x:importlib.metadata.version(x) for x in ('numberdb','mpmath')},
'rows':evidence}
print('EVIDENCE_JSON',json.dumps(summary),flush=True)
if __name__=='__main__':
if os.environ.get('NUMBERDB_KEY_FROM_STDIN')=='1':
os.environ['NUMBERDB_API_KEY']=sys.stdin.read().strip()
main('--stored' in sys.argv or os.environ.get('NUMBERDB_STORED')=='1')