back to table · edit · history · where entries came from · files · download
4047 bytes, as of the version from 2026-10-27 01:12 (current). Recorded here, not run.
"""Check both Coulomb auxiliary solutions -- numberdb.org/T511.
Run with SageMath and numberdb installed:
sage -python check-normalization.py
Keep T511-document.json beside this file. It is an ordered API snapshot.
This reads no credentials and writes nothing to NumberDB. DLMF 33.2.8
computes H+ and H- independently in 600-bit complex-ball arithmetic;
DLMF 33.2.11 is checked against mpmath F and G at 140 decimal digits.
Every resulting G ball must fit in the stored decimal interval. These
numerical checks supplement, rather than prove, the cited identity.
"""
import json
from decimal import Decimal
from pathlib import Path
import numberdb.sage as numberdb
from sage.rings.complex_arb import ComplexBallField
from sage.rings.rational_field import QQ
from sage.env import SAGE_VERSION
import mpmath as mp
C = ComplexBallField(600)
i = C(0, 1)
mp.mp.dps = 140
def rational_decimal(s):
sign, digs, exp = Decimal(s).as_tuple()
n = int(''.join(map(str, digs))) * (-1 if sign else 1)
return QQ(n) * QQ(10)**int(exp)
def as_mp(q):
q = QQ(q)
return mp.mpf(int(q.numerator())) / int(q.denominator())
def midpoint(x):
return mp.mpf(str(x.mid()))
def h(ell, eta, rho, sign):
a = C(ell + 1) + sign*i*eta
# The exponential needs only the phase modulo 2*pi; Gamma/|Gamma|
# therefore supplies it without choosing a logarithm branch.
gamma = a.gamma()
phase = (-sign*i)**ell * gamma/abs(gamma)
phase *= (sign*i*(rho-eta*(2*rho).log())).exp()
z = -sign*2*i*rho
return phase * z**a * z.hypergeometric_U(a, C(2*ell+2))
def narrow(z):
for x in (z.real(), z.imag()):
assert x.is_finite()
assert QQ(x.rad()) < QQ(1)/10**160
# A nonzero known answer checks the phase/sign convention first.
hp = h(0,C(0),C(1),1)
hm = h(0,C(0),C(1),-1)
assert hp.real().overlaps(C(1).real().cos())
assert hp.imag().overlaps(C(1).real().sin())
assert hm.imag().overlaps(-C(1).real().sin())
assert not hm.imag().overlaps(C(1).real().sin())
doc = json.loads(Path(__file__).with_name('T511-document.json').read_text())
checked = []
max_residual = mp.mpf(0)
max_ulps = mp.mpf(0)
for ls, etas in doc['Numbers'].items():
ell = int(ls)
for es, rhos in etas.items():
for rs, entry in rhos.items():
value = entry if isinstance(entry,str) else entry['number']
eta, rho = C(QQ(es)), C(QQ(rs))
hp, hm = h(ell,eta,rho,1), h(ell,eta,rho,-1)
narrow(hp); narrow(hm)
assert hp.conjugate().overlaps(hm)
g, f = (hp+hm)/2, (hp-hm)/(2*i)
assert g.imag().contains_zero() and f.imag().contains_zero()
center = rational_decimal(value)
unit = QQ(10)**int(Decimal(value).as_tuple().exponent)
assert g.real() > C(center-unit).real()
assert g.real() < C(center+unit).real()
fm = mp.coulombf(ell,as_mp(es),as_mp(rs))
gm = mp.coulombg(ell,as_mp(es),as_mp(rs))
for actual, expected in ((g.real(),gm),(f.real(),fm)):
residual = abs(midpoint(actual)-expected)/max(1,abs(expected))
assert residual < mp.mpf('1e-125')
max_residual = max(max_residual,residual)
max_ulps = max(max_ulps, abs(gm-as_mp(center))/as_mp(unit))
checked.append([ls,es,rs])
print('Checked order',ell,flush=True)
assert len(checked)==88
print(json.dumps({'table':'T511','entries':len(checked),'working_bits':600,
'mpmath_dps':mp.mp.dps,'sage_version':SAGE_VERSION,
'numberdb_version':numberdb.__version__,'mpmath_version':mp.__version__,
'max_scaled_mpmath_discrepancy':mp.nstr(max_residual,8),
'max_mpmath_discrepancy_in_last_place_units':mp.nstr(max_ulps,8),
'all_G_balls_inside_stored_intervals':True,'conjugacy_checked':True,
'free_case_sign_controls_passed':True,'reversed_sign_control_rejected':True,
'identities_source':'https://dlmf.nist.gov/33.2#E11',
'representation_source':'https://dlmf.nist.gov/33.2#E8',
'checked_parameters':checked},indent=2))