back to table · edit · history · where entries came from · files · download
6796 bytes, as of the version from 2026-10-27 01:12 (current). Recorded here, not run.
"""Zeros of the irregular Coulomb wave function -- numberdb.org/T513.
Run with SageMath or passagemath and numberdb:
sage -pip install numberdb mpmath
sage -python generate.py # verify all stored entries
sage -python checks.py # independent numerical checks
The script never publishes a table or submits it for review.
Convention: DLMF G_l(eta,rho), H+ = G+iF, rho>0, fixed rational eta.
The first strictly positive radial zero has index 1. Parameter order is
(ell, eta, n). The initial selection is ell=0..3, eta=0,+/-1/2,+/-1,+/-2,
and n=1..12, with the free case and smaller absolute couplings first.
Index certificate (restricted to |eta|<=2 and ell>=0):
Write F=C*rho**(ell+1)*sum(a_k*rho**k), where C>0,
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|<=5**k: the k=1 bound is 2<=5; for k>=2,
(4*5**(k-1)+5**(k-2))/(k*(k+2*ell+1))<=5**k.
Thus |sum_(k>=1)(a_k*rho**k)|<=5/11<1 for 0<rho<=1/16,
so F>0 there. By GF'-FG'=1, (G/F)'=-1/F**2<0.
A certified positive G(1/16) therefore proves G>0 throughout (0,1/16].
For rho>=1/16, q=1-2*eta/rho-ell*(ell+1)/rho**2<=65.
Between any two zeros a<b, integration by parts gives
integral(w'**2)=integral(q*w**2)<=65*integral(w**2).
The Dirichlet Wirtinger inequality gives
integral(w'**2)>=(pi/(b-a))**2*integral(w**2).
Consequently b-a>=pi/sqrt(65)>1/4. Zeros are simple by uniqueness
for the ODE. A scan in intervals of length 1/4 with strictly certified
endpoint signs therefore counts every zero, including the first.
The G evaluation uses DLMF 33.2.8 and principal U, with the powers
simplified for integer ell. Both H+ and H- are evaluated, and the imaginary
ball of their average must contain zero. At 160 bits, the scan checks
nonzero signs, increasing precision if necessary. Newton steps are only
candidate proposals. An answer is returned only after exact rational
endpoints, separated by 2*10**(-(digits+20)), are certified to have
opposite signs inside the corresponding indexed scan interval. Final
arithmetic uses numberdb.bits(digits+100), or 681 bits for 100 digits.
The returned interval encloses those rational endpoints, rather than a
fixed-precision number wrapped in an interval. checks.py also tests the
actual decimal spelling's endpoints and two-precision mpmath agreement.
Sources: https://dlmf.nist.gov/33.2, /33.4, /33.5 and /33.6.
"""
import os
import sys
from functools import lru_cache
import numberdb.sage as numberdb
from sage.rings.rational_field import QQ
from sage.rings.complex_arb import ComplexBallField
from sage.rings.real_mpfi import RealIntervalField
ETAS = ('0', '1/2', '-1/2', '1', '-1', '2', '-2')
START = QQ(1) / 16
STEP = QQ(1) / 4
SCAN_BITS = 160
def finite(z):
assert z.real().is_finite() and z.imag().is_finite(), 'nonfinite ball'
return z
@lru_cache(maxsize=None)
def prefactor(ell, eta, bits, sign):
C = ComplexBallField(bits)
i = C(0, sign)
e = C(QQ(eta))
a = C(ell + 1) + i*e
gamma = a.gamma()
# (-i)^ell from the phase AND (-i)^(ell+1) from the power.
pref = (-1)**ell * (-2*i) * C(2)**ell
pref *= (C.pi()*e/2).exp() * gamma/abs(gamma)
return finite(pref), a
def coulomb_g(ell, eta, rho, bits=681):
C = ComplexBallField(bits)
r = C(rho)
waves = []
for sign in (1, -1):
i = C(0, sign)
pref, a = prefactor(ell, str(eta), bits, sign)
h = pref*r**(ell+1)*(i*r).exp()
h *= (-2*i*r).hypergeometric_U(a, 2*ell+2)
waves.append(finite(h))
g = finite((waves[0]+waves[1])/2)
assert g.imag().contains_zero(), 'G must be real'
return g.real()
def sign_at(ell, eta, rho, bits=SCAN_BITS):
for prec in (bits, bits*2, bits*4):
v = coulomb_g(ell, eta, rho, prec)
if v > 0:
return 1
if v < 0:
return -1
raise ArithmeticError('could not certify nonzero endpoint sign')
@lru_cache(maxsize=None)
def scan(ell, eta, count=12):
assert ell >= 0 and abs(QQ(eta)) <= 2 and count >= 1
left = START
previous = sign_at(ell, eta, left)
assert previous == 1, 'initial interval proof requires G(1/16)>0'
brackets = []
for _ in range(10000):
right = left + STEP
current = sign_at(ell, eta, right)
if current != previous:
brackets.append((left, right))
if len(brackets) == count:
return tuple(brackets)
left, previous = right, current
raise ArithmeticError('scan limit reached')
@lru_cache(maxsize=None)
def certified_root(ell, eta, n, digits):
outer_a, outer_b = scan(ell, eta, max(12, n))[n-1]
a, b = outer_a, outer_b
sa = sign_at(ell, eta, a)
for _ in range(20):
m = (a+b)/2
if sign_at(ell, eta, m) == sa:
a = m
else:
b = m
q = (a+b)/2
bits = numberdb.bits(digits+100)
C = ComplexBallField(bits)
e = C(QQ(eta))
for _ in range(10):
g = coulomb_g(ell, eta, q, bits)
gn = coulomb_g(ell+1, eta, q, bits)
derivative = ((C(ell+1)/C(q)+e/(ell+1))*g
-(1+(e/(ell+1))**2).sqrt()*gn).real()
assert not derivative.contains_zero(), 'uncertain Newton derivative'
correction = g/derivative
q -= correction.mid().exact_rational()
assert outer_a < q < outer_b, 'Newton left indexed interval'
if abs(correction) < C(QQ(1)/10**(digits+30)).real():
break
radius = QQ(1)/10**(digits+20)
lower, upper = q-radius, q+radius
assert outer_a < lower < upper < outer_b
assert sign_at(ell, eta, lower, bits) == sa
assert sign_at(ell, eta, upper, bits) == -sa
return RealIntervalField(bits)(lower, upper)
class CoulombGZeros(numberdb.Generator):
table = 'T513'
parameters = ('ell', 'eta', 'n')
type = 'R'
digits = 100
rigour = 'proven'
def enumerate(self, max_ell=3, count=12):
for ell in range(max_ell+1):
for eta in ETAS:
for n in range(1, count+1):
yield {'ell': str(ell), 'eta': eta, 'n': str(n)}
def value(self, params, digits):
ell, eta, n = int(params['ell']), str(params['eta']), int(params['n'])
v = certified_root(ell, eta, n, digits)
result = {'number': v}
if QQ(eta) == 0:
order = str(QQ(2*ell+1)/2)
result['equals'] = ('HREF{Zeros_of_Bessel_Y_functions#%s,%d}'
'[Bessel zero $y_{%s,%d}$]' % (order, n, order, n))
return result
if __name__ == '__main__':
if os.environ.get('NUMBERDB_KEY_FROM_STDIN') == '1':
os.environ['NUMBERDB_API_KEY'] = sys.stdin.read().strip()
report = CoulombGZeros().verify(sample=None)
print(report)
sys.exit(0 if report.ok else 1)