back to table · edit · history · where entries came from · files · download
6650 bytes, as of the version from 2026-10-27 01:12 (current). Recorded here, not run.
"""Zeros of the regular Coulomb wave function -- numberdb.org/T512.
Run with SageMath (install numberdb and mpmath once):
sage -pip install numberdb mpmath
sage -python generate.py # verify the stored draft/table
sage -python generate.py --fill # fill the existing draft, key required
sage -python checks.py # independent checks, no key required
The 336 radial roots use ell=0..3, eta=0,+/-1/2,+/-1,+/-2, n=1..12.
A rigorous scan locates every root. mpmath at 150 digits supplies candidates;
Arb signs at exact rational endpoints certify 120-decimal-place brackets.
The output is a true Sage interval containing the indexed root, not an
interval wrapped around an uncertified numerical answer.
INDEX PROOF (valid for ell>=0, |eta|<=2): Write F=C*r**(ell+1)*sum a_k*r**k.
DLMF 33.6 gives a_0=1, a_1=eta/(ell+1) and
k*(k+2*ell+1)*a_k=2*eta*a_{k-1}-a_{k-2}. Induction gives |a_k|<=3**k:
the first coefficient is <=2, and for k>=2 the inductive multiplier is
(4/3+1/9)/(k*(k+2*ell+1)) <=13/54<1. Thus for 0<r<=1/8 the sum is
at least 1-sum((3*r)**k,k>=1)>=2/5, so there is no positive zero there.
For r>=1/8 the equation coefficient q(r)=1-2*eta/r-ell*(ell+1)/r**2
is <=33. Between two zeros a<b, integration by parts and the Dirichlet
Wirtinger inequality give (pi/(b-a))**2<=33. Hence their distance is
at least pi/sqrt(33)>1/2. All zeros are simple by ODE uniqueness.
Scanning exact points 1/8+j/4 with certified, nonzero signs therefore
cannot miss any root. Every sign change contains exactly one.
The scan uses 160 working bits with explicit failure if a sign is unresolved.
The narrow-bracket certification uses numberdb.bits(200)=681 bits; only
finite, real balls with strict opposite signs are accepted. A pilot at 515 bits
failed to resolve a 1e-120 endpoint sign near radius 40 because the Kummer
representation loses accuracy through cancellation. The precision
is a starting choice, not an accuracy claim: each accepted bracket has exact
rational width 2e-120 and is contained in its indexed scan interval.
See checks.py for controls, identities, independent precision agreement,
Bessel zeros, and certification of the decimals actually written.
"""
import os
import sys
from functools import lru_cache
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
import mpmath as mp
ETAS = ('0', '1/2', '-1/2', '1', '-1', '2', '-2')
WORK_BITS = numberdb.bits(200)
SCAN_BITS = 160
ENCLOSURE_PLACES = 120
def rational_decimal(s):
"""An exact rational from decimal digits, including scientific notation."""
s = str(s).strip().rstrip('?').lower()
mantissa, _, exponent = s.partition('e')
power = int(exponent or 0)
if '.' in mantissa:
whole, frac = mantissa.split('.')
power -= len(frac)
mantissa = whole + frac
return QQ(int(mantissa)) * QQ(10)**power
def written_interval(text):
s = str(text).strip().rstrip('?').lower()
m, _, e = s.partition('e')
ulp = QQ(10)**(int(e or 0) - (len(m.split('.')[1]) if '.' in m else 0))
center = rational_decimal(s)
return center-ulp, center+ulp
def regular_scaled(ell, eta, rho, bits=WORK_BITS):
"""F/(C*r**(ell+1)); its sign equals F's for positive real rho."""
C = ComplexBallField(bits)
i = C(0, 1)
r, e = C(rho), C(eta)
v = (-i*r).exp() * (2*i*r).hypergeometric([ell+1-i*e], [C(2*ell+2)])
assert v.real().is_finite() and v.imag().is_finite()
assert v.imag().contains_zero()
return v.real()
def sign_at(ell, eta, rho, bits=WORK_BITS):
v = regular_scaled(ell, eta, rho, bits)
if v > 0:
return 1
if v < 0:
return -1
raise ArithmeticError('Unresolved sign at %s,%s,%s; increase working precision' % (ell, eta, rho))
@lru_cache(None)
def scan(ell, eta, count=12):
ell, eta = int(ell), QQ(eta)
assert ell >= 0 and abs(eta) <= 2
a = QQ(1)/8
sa = sign_at(ell, eta, a, SCAN_BITS)
assert sa == 1
result = []
for j in range(1, 2001):
b = QQ(1)/8 + QQ(j)/4
sb = sign_at(ell, eta, b, SCAN_BITS)
if sa != sb:
result.append((a, b))
if len(result) == count:
return tuple(result)
a, sa = b, sb
raise ArithmeticError('Scan did not find the requested zeros')
def mpq(x):
x = QQ(x)
return mp.mpf(str(x.numerator())) / mp.mpf(str(x.denominator()))
@lru_cache(None)
def certified_root(ell, eta, n):
ell, n, eta = int(ell), int(n), QQ(eta)
assert 1 <= n <= 12
a, b = scan(ell, str(eta))[n-1]
with mp.workdps(150):
e = mpq(eta)
root = mp.findroot(lambda x: mp.coulombf(ell, e, x),
(mpq(a), mpq(b)), solver='secant', tol=mp.mpf('1e-140'))
candidate = rational_decimal(mp.nstr(root, 140))
radius = QQ(1)/10**ENCLOSURE_PLACES
lo, hi = candidate-radius, candidate+radius
assert a < lo < hi < b
assert sign_at(ell, eta, lo) == (-1)**(n-1)
assert sign_at(ell, eta, hi) == (-1)**n
return RealIntervalField(WORK_BITS)(lo, hi)
class CoulombFZeros(numberdb.Generator):
table = 'T512'
parameters = ('ell', 'eta', 'n')
type = 'R'
digits = 100
rigour = 'proven'
files = ('generate.py', 'checks.py', 'preflight-evidence.json')
def enumerate(self):
for ell in range(4):
for eta in ETAS:
for n in range(1, 13):
yield {'ell': str(ell), 'eta': eta, 'n': str(n)}
def value(self, params, digits):
if digits != 100:
raise ValueError('This certified run is designed for 100 significant digits')
ell, eta, n = int(params['ell']), QQ(params['eta']), int(params['n'])
value = certified_root(ell, str(eta), n)
if eta == 0:
return {'number': value,
'equals': ('HREF{Zeros_of_Bessel_J_functions#%s,%s}'
'[Bessel zero $j_{%s,%s}$]')
% (QQ(2*ell+1)/2, n, QQ(2*ell+1)/2, n)}
return value
if __name__ == '__main__':
if os.environ.get('NUMBERDB_KEY_FROM_STDIN') == '1':
os.environ['NUMBERDB_API_KEY'] = sys.stdin.read().strip()
g = CoulombFZeros()
if '--fill' in sys.argv or os.environ.get('NUMBERDB_PUBLISH') == '1':
print(g.publish(message='Enclose and count 336 positive Coulomb radial zeros', assisted_by='codex'))
else:
report = g.verify(sample=None)
print(report)
sys.exit(0 if report.ok else 1)