back to table · edit · history · where entries came from · files · download
10811 bytes, as of the version from 2026-10-27 01:12 (current). Recorded here, not run.
"""Independent checks for Coulomb G radial zeros -- numberdb.org/T513.
sage -python checks.py > checks.log
NUMBERDB_CHECK_STORED=1 sage -python checks.py > stored-checks.log
The second command reads the website document and repeats all ball identities,
actual decimal-endpoint signs, Bessel checks and interlacing. Reading a draft
requires NUMBERDB_API_KEY. The first run additionally compares independently
computed mpmath/Taylor Newton estimates at 150 and 200 decimal digits. The mpmath Coulomb functions are evaluated 1/8 beyond each root; an ODE
Taylor expansion transfers their values to the root neighborhood. This avoids
mpmath hypercomb convergence failure arbitrarily close to a free-case zero.
Estimates are a numerical check, not the source of the generator's proof.
The checking G implementation uses the unsimplified DLMF phase/power formula,
analytic U derivatives and the gamma phase. F uses the product normalization
and shifted M parameters. No finite differences or order recurrences produce
the derivatives that are used to test the order and derivative recurrences.
Each ENTRY line records the value actually checked. SUMMARY records counts,
precision bounds and versions. Controls run before accepting any check.
"""
import os, sys, json, time, platform
from math import factorial
from fractions import Fraction
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_mpfi import RealIntervalField
from sage.env import SAGE_VERSION
import mpmath as mp
from generate import CoulombGZeros, scan, coulomb_g
from check import exactness, measure, names_its_rings, agrees_with
BITS = numberdb.bits(200)
C = ComplexBallField(BITS)
I = C(0,1)
MAX_RADIUS = QQ(0)
def checked_real(z):
global MAX_RADIUS
assert z.real().is_finite() and z.imag().is_finite()
assert z.imag().contains_zero()
for part in (z.real(),z.imag()):
rad = part.rad().exact_rational()
assert rad < QQ(1)/10**130, 'wide ball'
MAX_RADIUS = max(MAX_RADIUS, rad)
return z.real()
def zero(z):
z = C(z)
assert checked_real(z).contains_zero(), 'identity failed'
def regular(ell, eta, rho):
r,e=C(rho),C(eta)
c2 = C(1) if eta == 0 else 2*C.pi()*e/((2*C.pi()*e).exp()-1)
for k in range(1,ell+1):
c2 *= e*e+k*k
norm=C(2)**ell*c2.sqrt()/factorial(2*ell+1)
a,b=C(ell+1)-I*e,C(2*ell+2)
z=2*I*r
ms=[z.hypergeometric([a+j],[b+j]) for j in range(3)]
h=norm*r**(ell+1)*(-I*r).exp()
d=C(ell+1)/r-I
m1=2*I*a/b*ms[1]
m2=-4*a*(a+1)/(b*(b+1))*ms[2]
return tuple(checked_real(v) for v in
(h*ms[0],h*(d*ms[0]+m1),
h*((d*d-C(ell+1)/(r*r))*ms[0]+2*d*m1+m2)))
def irregular(ell, eta, rho):
r,e=C(rho),C(eta)
a,b=C(ell+1)+I*e,C(2*ell+2)
sigma=a.log_gamma().imag()
theta=r-e*(2*r).log()-ell*C.pi()/2+sigma
z=-2*I*r
pref=(I*theta).exp()*z**a
us=[z.hypergeometric_U(a+j,b+j) for j in range(3)]
d=C(ell+1)/r+I
u1=2*I*a*us[1]
u2=-4*a*(a+1)*us[2]
hs=[pref*us[0],pref*(d*us[0]+u1),
pref*((d*d-C(ell+1)/(r*r))*us[0]+2*d*u1+u2)]
# Real part is G by definition. Check finite, narrow real enclosures.
out=[]
for h in hs:
assert h.real().is_finite() and h.imag().is_finite()
out.append(checked_real((h+h.conjugate())/2))
return tuple(out)
def decimal_interval(s):
f=Fraction(s)
mant,_,exp=s.lower().partition('e')
places=len(mant.split('.')[1]) if '.' in mant else 0
ulp=Fraction(10)**(int(exp or '0')-places)
return QQ(f.numerator)/f.denominator, QQ(ulp.numerator)/ulp.denominator
def strict_sign(v):
assert v.is_finite()
if v>0:return 1
if v<0:return -1
raise AssertionError('indeterminate sign')
def mpq(q):
return mp.mpf(int(q.numerator()))/int(q.denominator())
def mp_root_estimate(ell, eta, r, dps):
"""Independent G evaluation away from cancellation, then ODE Taylor solve."""
x,e=mpq(r+QQ(1)/8),mpq(eta)
g=mp.coulombg(ell,e,x)
gn=mp.coulombg(ell+1,e,x)
gp=((ell+1)/x+e/(ell+1))*g-mp.sqrt(1+(e/(ell+1))**2)*gn
a=[g,gp]
q0=x*x-2*e*x-ell*(ell+1)
q1=2*(x-e)
# From (x+t)^2 y'' + ((x+t)^2-2*eta*(x+t)-ell*(ell+1))*y=0.
# Counts differ between precision runs, checking truncation as well as rounding.
for k in range(8*dps):
numerator=2*x*k*(k+1)*a[k+1]+(k*(k-1)+q0)*a[k]
if k>=1:numerator+=q1*a[k-1]
if k>=2:numerator+=a[k-2]
a.append(-numerator/(x*x*(k+2)*(k+1)))
t=-mp.mpf(1)/8
for _ in range(4):
value=a[-1];derivative=mp.mpf(0)
for ak in reversed(a[:-1]):
derivative=derivative*t+value
value=value*t+ak
t-=value/derivative
return x+t,g,gp,x
def identities(ell, eta, r):
f,fp,fpp=regular(ell,eta,r)
g,gp,gpp=irregular(ell,eta,r)
q=C(1)-2*C(eta)/r-C(ell*(ell+1))/(r*r)
zero(fpp+q*f);zero(gpp+q*g)
zero(g*fp-f*gp-1)
zero(g-coulomb_g(ell,str(eta),r,BITS))
fn=regular(ell+1,eta,r)[0]
gn=irregular(ell+1,eta,r)[0]
rn=(1+(C(eta)/(ell+1))**2).sqrt()
sn=C(ell+1)/r+C(eta)/(ell+1)
zero(fp-sn*f+rn*fn);zero(gp-sn*g+rn*gn)
if ell:
fm=regular(ell-1,eta,r)[0];gm=irregular(ell-1,eta,r)[0]
rl=(1+(C(eta)/ell)**2).sqrt();sl=C(ell)/r+C(eta)/ell
zero(fp-rl*fm+sl*f);zero(gp-rl*gm+sl*g)
zero(rl*fm-(sl+sn)*f+rn*fn)
zero(rl*gm-(sl+sn)*g+rn*gn)
return f,fp,g,gp
def controls():
# A known nonzero trigonometric pair and Wronskian, plus rejection control.
f,fp,g,gp=identities(0,QQ(0),QQ(1))
zero(f-C(1).sin());zero(g-C(1).cos())
assert f>0 and g>0 and not (C(-g*fp+f*gp-1).real().contains_zero())
assert decimal_interval('3.14')==(QQ(314)/100,QQ(1)/100)
assert decimal_interval('12e2')==(QQ(1200),QQ(100))
# mpmath documentation, G and -G' at eta=2,rho=5.
published=[(1,'1.08148276','0.6028279961'),
(2,'1.496877075','0.5661803178'),
(3,'2.048694714','0.7959909551')]
for ell,gs,ds in published:
g,gp,_=irregular(ell,QQ(2),QQ(5))
for v,s in [(g,gs),(-gp,ds)]:
c,u=decimal_interval(s);assert abs(v-C(c).real())<C(u).real()
# Exact inequality underlying the initial series bound, with base cases.
for ell in range(4):
for eta in map(QQ,('0','1/2','-1/2','1','-1','2','-2')):
coeff=[QQ(1),eta/(ell+1)]
for k in range(2,30):
coeff.append((2*eta*coeff[-1]-coeff[-2])/(k*(k+2*ell+1)))
assert all(abs(a)<=5**k for k,a in enumerate(coeff))
assert (C.pi()/C(65).sqrt()).real()>C(QQ(1)/4).real()
with mp.workdps(150):
estimate,_,_,_=mp_root_estimate(0,QQ(0),QQ(157)/100,150)
assert abs(estimate-mp.pi/2)<mp.mpf('1e-140')
print('CONTROLS passed: nonzero sine/cosine, Wronskian, wrong-sign rejection, decimal intervals, six published values, coefficient bound',flush=True)
def main():
start=time.monotonic()
stored_mode=os.environ.get('NUMBERDB_CHECK_STORED')=='1'
if os.environ.get('NUMBERDB_KEY_FROM_STDIN')=='1':
os.environ['NUMBERDB_API_KEY']=sys.stdin.read().strip()
stored_mode=True
gen=CoulombGZeros()
doc=numberdb.table('T513') if stored_mode else None
controls()
assert not names_its_rings(__file__)
checked={};bessel={};mp_estimates={};nentries=0
regular_doc=numberdb.table('T512') if stored_mode else None
for params in gen.enumerate():
ell,eta,n=int(params['ell']),QQ(params['eta']),int(params['n'])
if stored_mode:
entry=doc['Numbers'][str(ell)][str(eta)][str(n)]
text=entry['number'] if isinstance(entry,dict) else entry
else:
value=gen.value(params,100)
text=to_text(value['number'],100)
entry={**value,'number':text}
r,u=decimal_interval(text);lo,hi=r-u,r+u
f,fp,g,gp=identities(ell,eta,r)
ga=irregular(ell,eta,lo)[0];gb=irregular(ell,eta,hi)[0]
sa,sb=strict_sign(ga),strict_sign(gb)
assert sa==-sb
a,b=scan(ell,str(eta))[n-1]
assert a<lo<hi<b
# Scan starts positive and each successive simple zero reverses sign.
assert sa==(-1)**(n-1)
key=(ell,str(eta),n)
checked[key]=text
if eta==0:
with mp.workdps(200):
br=mp.besselyzero(mp.mpf(ell)+mp.mpf('0.5'),n)
assert mpq(lo)<br<mpq(hi)
if ell==0:assert abs(br-(mp.mpf(n)-mp.mpf('0.5'))*mp.pi)<mp.mpf('1e-195')
bessel[key]=text
record={'params':params,'number':text,'scan':[str(a),str(b)],'signs':[sa,sb]}
if not stored_mode:
estimates=[]
for dps in (150,200):
with mp.workdps(dps):
estimate,mg,dg,x=mp_root_estimate(ell,eta,r,dps)
assert mpq(lo)<estimate<mpq(hi)
estimates.append(mp.nstr(estimate,145))
ag,ad,_=irregular(ell,eta,r+QQ(1)/8)
assert abs(mg-mpq(ag.mid().exact_rational()))<mpq(ag.rad().exact_rational())+mp.mpf(10)**(-dps+8)
assert abs(dg-mpq(ad.mid().exact_rational()))<mpq(ad.rad().exact_rational())+mp.mpf(10)**(-dps+8)
with mp.workdps(200):
assert abs(mp.mpf(estimates[0])-mp.mpf(estimates[1]))<mp.mpf('1e-140')
record['mpmath_newton_150_200']=estimates
else:
# From separately stored F zeros: G_n < F_n < G_(n+1).
fe=regular_doc['Numbers'][str(ell)][str(eta)][str(n)]
fr,fu=decimal_interval(fe['number'] if isinstance(fe,dict) else fe)
assert hi<fr-fu
if n>1:
prev=regular_doc['Numbers'][str(ell)][str(eta)][str(n-1)]
pr,pu=decimal_interval(prev['number'] if isinstance(prev,dict) else prev)
assert pr+pu<lo
nentries+=1
print('ENTRY '+json.dumps(record),flush=True)
if n==12:print('PROGRESS',ell,str(eta),'entries',nentries,'seconds',round(time.monotonic()-start,1),flush=True)
assert nentries==336
assert not exactness(checked)
summary={'table':'T513','stored':stored_mode,'entries':nentries,'bessel_entries':len(bessel),
'bits':BITS,'max_function_ball_radius':str(C(MAX_RADIUS).real()),
'size':measure(checked),'mpmath':mp.__version__,'sage':SAGE_VERSION,
'python':platform.python_version(),'seconds':round(time.monotonic()-start,2),
'checks':'all decimal endpoint brackets, indices, ODE, Wronskian, analytic derivatives and order recurrences; Bessel and cosine controls; interlacing on stored run'}
print('SUMMARY '+json.dumps(summary),flush=True)
if __name__=='__main__':main()