back to table · edit · history · where entries came from · files · download
9939 bytes, as of the version from 2026-10-26 09:23. Recorded here, not run.
"""Independent checks for regular Coulomb values -- numberdb.org/T510.
Run sage -python checks.py for local checks. Run with --stored (and a key
allowed to read T510) to check all read-back decimal intervals. In the campaign,
NUMBERDB_KEY_FROM_STDIN=1 also selects stored mode and consumes the key on stdin.
The generator evaluates the M representation. This checker differentiates M
analytically with a separate product normalization, evaluates H+ through U,
and compares F with mpmath at 150 and 200 decimal digits. Bessel J and sine
supply checks outside the Coulomb family. Printed A&S-example values from the
mpmath documentation supply four low-precision F and F' controls. The original
A&S p.543 Example 1 is checked separately to absolute tolerance 1e-4: its
printed five-digit normalization loses accuracy in the final displayed digits.
Ball
residuals must contain zero AND have bounded absolute size; a wide or NaN ball
cannot pass. Stored decimals are parsed as exact endpoints at +/- one ulp.
"""
import json
import os
import sys
from decimal import Decimal
from functools import lru_cache
from math import factorial
from pathlib import Path
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.env import SAGE_VERSION
import mpmath as mp
from generate import RegularCoulomb
from check import exactness, measure, names_its_rings
from numberdb._write import to_text
from numberdb._generate import _stored_entries
C = ComplexBallField(600)
R = RealBallField(600)
I = C(0, 1)
def exact_decimal(text):
n, d = Decimal(text).as_integer_ratio()
return QQ(n) / QQ(d)
def decimal_interval(text):
d = Decimal(text)
centre = exact_decimal(text)
ulp = QQ(10)**int(d.as_tuple().exponent)
return centre - ulp, centre + ulp
def interval_ball(text):
lo, hi = decimal_interval(text)
return R(RealIntervalField(600)(lo, hi))
def finite_narrow(value, tolerance=QQ(1)/10**120):
assert value.is_finite(), 'nonfinite ball'
assert value.rad() < tolerance, 'wide ball'
def zero(value, tolerance):
finite_narrow(value, tolerance)
assert value.contains_zero(), 'identity residual excludes zero'
assert abs(value).upper() < tolerance, 'identity residual too large'
def real_part(value):
assert value.real().is_finite() and value.imag().is_finite()
zero(value.imag(), QQ(1)/10**150)
return value.real()
def mp_exact(value):
sign, mantissa, exponent, _ = value._mpf_
return (-1 if sign else 1) * QQ(mantissa) * QQ(2)**exponent
def mp_rational(value):
q = QQ(value)
return mp.mpf(int(q.numerator())) / int(q.denominator())
@lru_cache(None)
def auxiliary(ell, eta, rho):
"""Independent F derivatives, G and G' (no generator helper calls)."""
ell = int(ell)
e, r = R(QQ(eta)), R(QQ(rho))
if QQ(eta) == 0:
normal = R(QQ(2)**ell * factorial(ell) / factorial(2*ell+1))
else:
normal2 = 2*R.pi()*e / (2*R.pi()*e).expm1()
for k in range(1, ell+1):
normal2 *= e*e + k*k
normal = 2**ell * normal2.sqrt() / factorial(2*ell+1)
a, b, z = ell+1-I*e, C(2*ell+2), 2*I*r
m0 = z.hypergeometric([a], [b])
dm = 2*I*a/b * z.hypergeometric([a+1], [b+1])
ddm = -4*a*(a+1)/(b*(b+1)) * z.hypergeometric([a+2], [b+2])
pref = normal * r**(ell+1) * (-I*r).exp()
slope = (ell+1)/r-I
f = real_part(pref*m0)
fp = real_part(pref*(slope*m0+dm))
fpp = real_part(pref*((slope*slope-(ell+1)/(r*r))*m0+2*slope*dm+ddm))
# DLMF 33.2.7 with W expressed through U. Gamma/abs(Gamma) is
# exp(i sigma) without making a numerical phase-unwrapping decision.
aplus, zminus = ell+1+I*e, -2*I*r
gamma = aplus.gamma()
hpref = (-I)**ell * (R.pi()*e/2).exp() * gamma/abs(gamma) \
* (I*r).exp() * zminus**(ell+1)
u = zminus.hypergeometric_U(aplus, b)
du = 2*I*aplus * zminus.hypergeometric_U(aplus+1, b+1)
h = hpref*u
hp = hpref*((I+(ell+1)/r)*u+du)
for value in [f, fp, fpp, h.real(), h.imag(), hp.real(), hp.imag()]:
finite_narrow(value, QQ(1)/10**150)
zero(h.imag()-f, QQ(1)/10**150)
zero(hp.imag()-fp, QQ(1)/10**150)
return f, fp, fpp, h.real(), hp.real()
def main():
keyed = os.environ.get('NUMBERDB_KEY_FROM_STDIN') == '1'
if keyed:
os.environ['NUMBERDB_API_KEY'] = sys.stdin.read().strip()
stored_mode = keyed or '--stored' in sys.argv
generator = RegularCoulomb()
assert not names_its_rings(str(Path(__file__).with_name('generate.py')))
params = list(generator.enumerate())
balls = {','.join(p[k] for k in generator.parameters): generator.value(p,100) for p in params}
assert len(balls) == 88 and not exactness(balls)
written = {key: to_text(value,100) for key,value in balls.items()}
if stored_mode:
document = numberdb.table('T510')
written = _stored_entries(document)
assert set(written) == set(balls), 'stored parameter set differs'
# Control the comparison machinery before trusting any agreement:
lo, hi = decimal_interval('1.0000000000')
assert lo == QQ(9999999999)/10**10 and hi == QQ(10000000001)/10**10
assert not (lo < QQ(100000001)/10**8 < hi), 'incorrect value passed control'
known_f, known_fp, _, known_g, known_gp = auxiliary(0,'0','1')
zero(known_f-R(1).sin(),QQ(1)/10**150)
zero(known_g-R(1).cos(),QQ(1)/10**150)
known_w = known_g*known_fp-known_f*known_gp
zero(known_w-1,QQ(1)/10**150)
assert not (known_w-2).contains_zero(), 'Wronskian control is vacuous'
print('CONTROLS: sine, cosine, Wronskian=1, wrong-value rejection passed',flush=True)
tol = QQ(1)/10**85 if stored_mode else QQ(1)/10**120
used = {key: interval_ball(text) for key,text in written.items()} if stored_mode else balls
counts = {'mpmath_150':0,'mpmath_200':0,'ode':0,'wronskian':0,'upward_derivative':0,
'downward_derivative':0,'order_recurrence':0,'bessel':0,'sine':0,'printed_control':0,'handbook_control':0}
printed_f = ['1.143337392','0.9935056752','0.7212774133','0.4313159311']
printed_fp = ['0.2937960375','0.3708676452','0.3952408216','0.320534053']
handbook_f = ['1.1433','0.99346','0.72124','0.43130']
handbook_fp = ['0.29380','0.3709','0.3952','0.3205']
handbook_scaled = ['6591.81','5727.97','4158.46','2486.72']
# A&S p.543, Example 1, column 2 uses the rounded normalization
# 1.7344e-4. Reproduce that arithmetic; it is not a five-digit oracle.
for index in range(1,4):
normalized = exact_decimal(handbook_scaled[index])*exact_decimal('1.7344e-4')
assert abs(normalized-exact_decimal(handbook_f[index])) < QQ(1)/10**5
for p in params:
ell,e,r = int(p['ell']),p['eta'],p['rho']
key=','.join(p[k] for k in generator.parameters)
f=used[key]
finite_narrow(f,tol)
lo,hi=decimal_interval(written[key])
# Compare exact rational decimal endpoints with outward dyadic bounds.
assert lo <= QQ(balls[key].lower()) and QQ(balls[key].upper()) <= hi
for dps in (150,200):
with mp.workdps(dps):
external=mp.coulombf(ell,mp_rational(e),mp_rational(r))
assert mp.isfinite(external) and not isinstance(external,mp.mpc)
external_q=mp_exact(external)
assert lo < external_q < hi, 'mpmath disagrees with written interval'
assert QQ(abs(balls[key]-R(external_q)).upper()) < abs(external_q)/10**140
counts['mpmath_'+str(dps)]+=1
f_ind, fp, fpp, g, gp=auxiliary(ell,e,r)
zero(f-f_ind,tol)
rho,eta=R(QQ(r)),R(QQ(e))
zero(fpp+(1-2*eta/rho-ell*(ell+1)/(rho*rho))*f,tol)
counts['ode']+=1
zero(g*fp-f*gp-1,tol)
counts['wronskian']+=1
n=ell+1
sn=n/rho+eta/n
rn=(1+(eta/n)**2).sqrt()
next_key=f'{n},{e},{r}'
next_f=used[next_key] if next_key in used else auxiliary(n,e,r)[0]
zero(fp-sn*f+rn*next_f,tol)
counts['upward_derivative']+=1
if ell>0:
sl=ell/rho+eta/ell
rl=(1+(eta/ell)**2).sqrt()
previous_f=used[f'{ell-1},{e},{r}']
zero(fp-rl*previous_f+sl*f,tol)
counts['downward_derivative']+=1
zero(rl*previous_f-(sl+sn)*f+rn*next_f,tol)
counts['order_recurrence']+=1
if QQ(e)==0:
bessel=real_part((C.pi()*rho/2).sqrt()*C(rho).bessel_J(C(QQ(2*ell+1)/2)))
zero(f-bessel,tol)
counts['bessel']+=1
if ell==0:
zero(f-rho.sin(),tol)
counts['sine']+=1
if e=='2' and r=='5':
for value,text in [(f,printed_f[ell]),(fp,printed_fp[ell])]:
plo,phi=decimal_interval(text)
assert plo < QQ(value.lower()) and QQ(value.upper()) < phi
counts['printed_control']+=1
for value,text in [(f,handbook_f[ell]),(fp,handbook_fp[ell])]:
assert QQ(abs(value-R(exact_decimal(text))).upper()) < QQ(1)/10**4
counts['handbook_control']+=1
print('ENTRY '+json.dumps({'params':p,'number':written[key],
'accuracy_bits':int(balls[key].accuracy()),'radius':str(balls[key].rad())}),flush=True)
summary={'mode':'stored' if stored_mode else 'local','tid':'T510',
'counts':counts,'entries':len(balls),'minimum_accuracy_bits':min(int(x.accuracy()) for x in balls.values()),
'maximum_radius':str(max(x.rad() for x in balls.values())),
'written_measure':measure(written),'Sage':SAGE_VERSION,'mpmath':mp.__version__,
'numberdb':getattr(numberdb,'__version__','unknown'),'passed':True}
print('SUMMARY '+json.dumps(summary,default=str),flush=True)
if __name__=='__main__':
main()