"""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 00 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=(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)