back to table · edit · history · where entries came from · files · download
10056 bytes, as of the version from 2026-10-26 10:01. Recorded here, not run.
"""Numerical checks for irregular Coulomb values -- numberdb.org/T511.
sage -python checks.py # all local checks, no credentials
sage -python checks.py --stored # check database values; draft needs key
Uses independent mpmath Coulomb G at 150 and 200 decimal digits; analytic
U derivatives; M-based regular solutions with a product normalization;
Bessel Y and cosine; and transcribed A&S / mpmath numerical controls.
The reported JSON is reproducible evidence, not a certificate from mpmath.
The value certification is from Arb enclosures and exact decimal endpoints.
"""
import numberdb.sage as numberdb
from sage.rings.complex_arb import ComplexBallField
from sage.rings.real_arb import RealBallField
from sage.rings.real_mpfi import RealIntervalField
from sage.rings.rational_field import QQ
from sage.version import version as sage_version
from numberdb._write import to_text
from numberdb._generate import _stored_entries
from generate import IrregularCoulomb
from decimal import Decimal
from fractions import Fraction
from math import factorial
import mpmath as mp
import json
import os
import sys
import time
BITS = numberdb.bits(180)
C = ComplexBallField(BITS)
R = RealBallField(BITS)
RI = RealIntervalField(BITS)
def exact_decimal(text):
return QQ(str(Fraction(Decimal(str(text).rstrip('?')))))
def interval_of(text):
d = Decimal(str(text).rstrip('?'))
center = QQ(str(Fraction(d)))
ulp = QQ(10)**int(d.as_tuple().exponent)
return RI(center-ulp, center+ulp), center, ulp
def bounded(x, scale=1, places=130):
assert x.is_finite(), 'nonfinite enclosure'
assert x.rad_as_ball() < R(10)**(-places)*max(R(1),abs(R(scale))), 'wide enclosure'
def zero(x, scale=1, places=130):
bounded(x,scale,places)
assert x.contains_zero(), 'identity residual excludes zero: %s' % x
def regular(ell, eta, rho):
"""Independent regular solution and derivatives, using a product for C_l."""
e,r,i=C(QQ(eta)),C(QQ(rho)),C(0,1)
if QQ(eta)==0:
norm=C(1)
else:
norm=(2*C.pi()*e/((2*C.pi()*e).exp()-1)).sqrt()
for j in range(1,ell+1):
norm *= (j*j+e*e).sqrt()
norm *= C(2)**ell / factorial(2*ell+1)
a,b,z=ell+1-i*e,C(2*ell+2),2*i*r
m=z.hypergeometric([a],[b])
m1=z.hypergeometric([a+1],[b+1])
m2=z.hypergeometric([a+2],[b+2])
pref=norm*r**(ell+1)*(-i*r).exp()
A=(ell+1)/r-i
f=pref*m
df=pref*(A*m+2*i*a/b*m1)
ddf=pref*((A*A-(ell+1)/(r*r))*m+4*i*A*a/b*m1-4*a*(a+1)/(b*(b+1))*m2)
for x in (f,df,ddf): zero(x.imag(),abs(x.real()))
return f.real(),df.real(),ddf.real()
def irregular_derivatives(ell,eta,rho):
"""U derivatives and ODE residual independent of order recurrences."""
e,r,i=C(QQ(eta)),C(QQ(rho)),C(0,1)
a,b,z=ell+1+i*e,C(2*ell+2),-2*i*r
gamma=a.gamma()
pref=(-i)**ell*(C.pi()*e/2).exp()*gamma/abs(gamma)*(i*r).exp()*z**(ell+1)
u=z.hypergeometric_U(a,b)
u1=z.hypergeometric_U(a+1,b+1)
u2=z.hypergeometric_U(a+2,b+2)
A=(ell+1)/r+i
h=pref*u
theta=r-e*(2*r).log()-ell*C.pi()/2+a.log_gamma().imag()
h_dlmf=(i*theta).exp()*z**a*u
zero((h-h_dlmf).real(),abs(h.real()))
zero((h-h_dlmf).imag(),abs(h.imag()))
dh=pref*(A*u+2*i*a*u1)
ddh=pref*((A*A-(ell+1)/(r*r))*u+4*i*A*a*u1-4*a*(a+1)*u2)
return h.real(),dh.real(),ddh.real(),h.imag(),dh.imag()
def run(stored=False):
started=time.monotonic()
gen=IrregularCoulomb()
if stored:
doc=numberdb.table(gen.table)
texts={tuple(identity.split(',')):text for identity,text in _stored_entries(doc).items()}
expected={tuple(p[k] for k in gen.parameters) for p in gen.enumerate()}
assert set(texts)==expected
else:
texts={}
# Known nonzero controls precede comparisons; wrong sign must fail.
g0=gen.value({'ell':'0','eta':'0','rho':'1'},100)
zero(g0-R(1).cos())
assert not (g0+R(1).cos()).contains_zero()
f0,df0,ddf0=regular(0,'0','1')
zero(f0-R(1).sin());zero(df0-R(1).cos());zero(ddf0+R(1).sin())
_,dg0,_,_,_=irregular_derivatives(0,'0','1')
w0=g0*df0-f0*dg0
zero(w0-1); assert not w0.contains_zero()
# Controls for the decimal parser: 1.25 denotes [1.24,1.26].
control,center,ulp=interval_of('1.25')
assert center==QQ(5)/4 and ulp==QQ(1)/100
assert QQ(124)/100 in control and QQ(123)/100 not in control
counts={'decimal_enclosures':0,'mpmath_150':0,'mpmath_200':0,'precision_agreement':0,
'ode_G':0,'ode_F':0,'wronskian':0,'order_G':0,'order_F':0,
'derivative_up_G':0,'derivative_down_G':0,'bessel_Y':0,'cosine':0,
'handbook_values':0,'handbook_derivatives':0,'mpmath_printed':0,
'Hplus_imag_F':0}
values={}; minbits=100000;maxrad=R(0); max_ulp_fraction=mp.mpf(0)
handbook_g=('0.79445','1.0815','1.4969','2.0487')
handbook_dg=('-0.67049','-0.60286','-0.56619','-0.79597')
printed_g={1:'1.08148276',2:'1.496877075',3:'2.048694714'}
printed_mdg={1:'0.6028279961',2:'0.5661803178',3:'0.7959909551'}
for index,p in enumerate(gen.enumerate(),1):
key=tuple(p[k] for k in gen.parameters)
ell=int(p['ell']);e=R(QQ(p['eta']));r=R(QQ(p['rho']))
local=gen.value(p,100)
bounded(local,local)
minbits=min(minbits,int(local.accuracy()));maxrad=max(maxrad,local.rad_as_ball())
written=texts[key] if stored else to_text(local,digits=100)
enclosure,center,ulp=interval_of(written)
local_interval=RI(local)
assert local_interval.lower()>=enclosure.lower() and local_interval.upper()<=enclosure.upper()
counts['decimal_enclosures']+=1
g=R(enclosure) if stored else local
scale=max(R(1),abs(g))
places=90 if stored else 125
values[key]=g
f,df,ddf=regular(ell,p['eta'],p['rho'])
gu,dg,ddg,imag_h,imag_dh=irregular_derivatives(ell,p['eta'],p['rho'])
zero(g-gu,scale,places)
zero(imag_h-f,max(abs(f),scale),125);zero(imag_dh-df,max(abs(df),scale),125)
counts['Hplus_imag_F']+=1
coeff=1-2*e/r-ell*(ell+1)/(r*r)
zero(ddg+coeff*g,scale*max(1,abs(coeff)),places);counts['ode_G']+=1
zero(ddf+coeff*f,max(1,abs(f)*abs(coeff)),125);counts['ode_F']+=1
zero(g*df-f*dg-1,max(1,abs(g*df),abs(f*dg)),places);counts['wronskian']+=1
def rr(j): return (1+e*e/(j*j)).sqrt()
def ss(j): return j/r+e/j
pup={'ell':str(ell+1),'eta':p['eta'],'rho':p['rho']}
upkey=tuple(pup[k] for k in gen.parameters)
if stored and upkey in texts:
gp=R(interval_of(texts[upkey])[0])
else: gp=gen.value(pup,100)
fp=regular(ell+1,p['eta'],p['rho'])[0]
zero(dg-ss(ell+1)*g+rr(ell+1)*gp,max(scale*abs(ss(ell+1)),abs(gp)),places)
counts['derivative_up_G']+=1
if ell:
downkey=(str(ell-1),p['eta'],p['rho'])
gm=values[downkey]
fm=regular(ell-1,p['eta'],p['rho'])[0]
zero(rr(ell)*gm-(ss(ell)+ss(ell+1))*g+rr(ell+1)*gp,max(abs(gm),abs(gp),scale*abs(ss(ell)+ss(ell+1))),places)
counts['order_G']+=1
zero(rr(ell)*fm-(ss(ell)+ss(ell+1))*f+rr(ell+1)*fp,max(1,abs(fm),abs(fp)),125)
counts['order_F']+=1
zero(dg-rr(ell)*gm+ss(ell)*g,max(scale*abs(ss(ell)),abs(gm)),places)
counts['derivative_down_G']+=1
if p['eta']=='0':
b=-(R.pi()*r/2).sqrt()*C(r).bessel_Y(C(QQ(ell)+QQ(1)/2)).real()
bounded(b,scale)
zero(g-b,scale,places);counts['bessel_Y']+=1
if ell==0:zero(g-r.cos(),scale,places);counts['cosine']+=1
if p['eta']=='2' and p['rho']=='5':
assert abs(gu-R(exact_decimal(handbook_g[ell])))<R(QQ(1)/10000)
assert abs(dg-R(exact_decimal(handbook_dg[ell])))<R(QQ(1)/10000)
counts['handbook_values']+=1;counts['handbook_derivatives']+=1
if ell in printed_g:
for actual,s in ((gu,printed_g[ell]),(-dg,printed_mdg[ell])):
target,c,u=interval_of(s)
assert abs(actual-R(c))<R(u)
counts['mpmath_printed']+=1
answers=[]
for dps in (150,200):
with mp.workdps(dps):
ee=mp.mpf(str(QQ(p['eta']).numerator()))/mp.mpf(str(QQ(p['eta']).denominator()))
rad=mp.mpf(str(QQ(p['rho']).numerator()))/mp.mpf(str(QQ(p['rho']).denominator()))
ans=mp.coulombg(ell,ee,rad)
assert mp.isfinite(ans)
mid=mp.mpf(str(center.numerator()))/mp.mpf(str(center.denominator()))
err=mp.mpf(str(ulp.numerator()))/mp.mpf(str(ulp.denominator()))
fraction=abs(ans-mid)/err
assert fraction<1
max_ulp_fraction=max(max_ulp_fraction,fraction)
answers.append(ans)
counts['mpmath_'+str(dps)]+=1
with mp.workdps(210):
assert abs(answers[0]-answers[1])<mp.mpf('1e-145')*max(1,abs(answers[1]))
counts['precision_agreement']+=1
if index%11==0: print('checked',index,'of 88',flush=True)
report={'table':gen.table,'mode':'stored' if stored else 'local','counts':counts,
'minimum_relative_accuracy_bits':minbits,'maximum_ball_radius':str(maxrad),
'working_bits':numberdb.bits(150),'check_bits':BITS,
'max_mpmath_distance_in_written_ulps':mp.nstr(max_ulp_fraction,8),
'entries':len(values),'versions':{'sage':sage_version,'numberdb':numberdb.__version__,'mpmath':mp.__version__},
'controls':['G0(0,1)=cos(1) and opposite sign rejected','F0(0,1)=sin(1)',
'Wronskian encloses 1 and excludes 0','decimal 1.25 parsed as [1.24,1.26]'],
'seconds':round(time.monotonic()-started,2)}
print('EVIDENCE_JSON='+json.dumps(report,sort_keys=True))
return report
if __name__=='__main__':
if os.environ.get('NUMBERDB_KEY_FROM_STDIN')=='1':
os.environ['NUMBERDB_API_KEY']=sys.stdin.read().strip()
run(stored='--stored' in sys.argv)