back to table · edit · history · where entries came from · files · download
9477 bytes, as of the version from 2026-10-25 13:21 (current). Recorded here, not run.
"""Zeros of cross-products of Bessel functions -- numberdb.org/T487
Run with SageMath (or a modular Sage installation providing these rings):
sage -pip install numberdb mpmath scipy
sage -python generate.py # recompute and verify the draft
sage -python generate.py --compute # local numerical evidence only
sage -python generate.py --fill # write entries to the draft
Authenticated access to this unpublished table needs NUMBERDB_API_KEY.
--fill sends numbers and source; it does not publish or offer the table.
The campaign wrapper can instead supply the key on stdin and set
NUMBERDB_KEY_FROM_STDIN=1 and NUMBERDB_PUBLISH=1.
Convention: k-th positive zero of J_nu(x)Y_nu(lambda*x)-Y_nu(x)J_nu(lambda*x).
The inner radius is 1. The derivative and mixed cross-products are excluded.
Certification, including the index:
* SciPy supplies guesses only. mpmath recomputes them at digits+30 and
digits+60 decimal places; agreement to digits+20 places is required.
* Exact rational endpoints at distance 10^(-digits-10) from a decimal
centre have strictly opposite signs certified by Arb. The returned ball
encloses both endpoints. No floating-point approximation is used as a bound.
* Ten ordered, disjoint sign brackets give at least ten roots. To rule out
missed roots, count radial nodes at the last bracket's upper endpoint X.
u(r)=sqrt(r)*(J_nu(X)Y_nu(r*X)-Y_nu(X)J_nu(r*X)) satisfies
-u''+(nu^2-1/4)*u/r^2=X^2*u, u(1)=0, u'(1)=2/pi>0.
Its consecutive zeros have distance at least pi/sqrt(X^2+1/4): integrate
the equation between two zeros and apply the Dirichlet Poincare inequality.
Choose an exact rational mesh h with h^2*(X^2+1/4)<9<pi^2.
Every mesh interval contains at most one zero; strict Arb sample signs
therefore count all interior nodes. The initial positive sign follows
from u'(1)>0 and the same spacing bound from the zero at r=1.
Sturm oscillation (Teschl, Theorem 5.18) identifies this node count with
the number of Dirichlet eigenvalues below X^2. All eigenvalues are positive,
by the original radial energy integral integral(r*R'^2+nu^2*R^2/r).
Exactly ten nodes proves there are exactly ten roots below X; hence each
bracket contains the correctly indexed, unique root.
Arb works at numberdb.bits(digits, losing=384) bits. At 100 digits this
is 717 bits; the enclosures are about 110 decimal places wide. The guard accounts for about 70 decimal digits lost by complex Arb Bessel
evaluation near the largest arguments (ratio=5/4, k=10). The measured
accuracy and every rational endpoint are emitted by --compute.
"""
import numberdb.sage as numberdb
from sage.rings.rational_field import QQ
from sage.rings.real_arb import RealBallField
from sage.rings.complex_arb import ComplexBallField
from functools import lru_cache
import json
import math
import os
import sys
import time
import mpmath as mp
from scipy.optimize import brentq
from scipy.special import jv, yv
RATIOS = ('5/4', '3/2', '2', '5/2', '3', '4', '5', '10')
MAX_ORDER = 5
ROOTS = 10
GUARD_BITS = 384
def decimal_rational(text):
"""Read a finite decimal exactly without Sage's optional decimal parser."""
mantissa, _, exponent = str(text).lower().partition('e')
scale = int(exponent or '0')
whole, dot, fraction = mantissa.partition('.')
return QQ(int(whole + fraction)) * QQ(10)**(scale-len(fraction))
def strict_sign(ball):
if not ball.is_finite():
raise ArithmeticError('non-finite ball')
if ball > 0:
return 1
if ball < 0:
return -1
raise ArithmeticError('ball does not certify a strict sign')
def cross_product(nu, ratio, x, bits):
C = ComplexBallField(bits)
x, ratio, nu = C(x), C(ratio), C(nu)
j, y = x.bessel_J_Y(nu)
jl, yl = (ratio*x).bessel_J_Y(nu)
value = j*yl-y*jl
if not value.real().is_finite() or not value.imag().is_finite() or not value.imag().contains_zero():
raise ArithmeticError('cross-product is not finite and real')
return value.real()
def radial_node_count(nu, ratio, x, bits):
"""Certified Sturm count for a positive rational x which is not a root."""
nu, ratio, x = QQ(nu), QQ(ratio), QQ(x)
if nu < 0 or ratio <= 1 or x <= 0:
raise ValueError('nu>=0, ratio>1 and x>0 are required')
m = int(((ratio-1)*(x+1)/QQ(2)).ceil()) + 1
h = (ratio-1)/m
assert h*h*(x*x+QQ(1)/4) < 9
C = ComplexBallField(bits)
j, y = C(x).bessel_J_Y(C(nu))
previous, count = 1, 0
signs = []
for i in range(1, m+1):
r = 1+i*h
jr, yr = C(r*x).bessel_J_Y(C(nu))
f = j*yr-y*jr
assert f.real().is_finite() and f.imag().is_finite() and f.imag().contains_zero()
sign = strict_sign(f.real())
if sign != previous:
count += 1
signs.append(sign)
previous = sign
return count, {'mesh_intervals':m, 'mesh_width':str(h), 'signs':signs}
def approximate_roots(nu, ratio, count):
"""Non-rigorous candidates; certification later detects missing roots."""
order, lam = float(nu), float(ratio)
def f(x):
return jv(order,x)*yv(order,lam*x)-yv(order,x)*jv(order,lam*x)
step = math.pi/(8*(lam-1))
a = min(1e-6,step/100)
fa = f(a)
found = []
for _ in range(10000):
b = a+step
fb = f(b)
if not math.isfinite(fa) or not math.isfinite(fb):
raise ArithmeticError('non-finite initial scan')
if fa*fb < 0:
found.append(brentq(f,a,b,xtol=1e-14))
if len(found) == count:
return found
a, fa = b, fb
raise ArithmeticError('initial scan exhausted')
def high_precision_root(nu, ratio, guess, digits):
with mp.workdps(digits):
v = mp.mpf(str(nu.numerator()))/mp.mpf(str(nu.denominator()))
lam = mp.mpf(str(ratio.numerator()))/mp.mpf(str(ratio.denominator()))
def f(x):
return mp.besselj(v,x)*mp.bessely(v,lam*x)-mp.bessely(v,x)*mp.besselj(v,lam*x)
g = mp.mpf(repr(guess))
shift = mp.mpf('0.00001')/(lam-1)
root = mp.findroot(f,(g-shift,g+shift),solver='secant',maxsteps=50)
if not mp.isfinite(root) or root <= 0:
raise ArithmeticError('invalid high-precision root')
return mp.nstr(root,digits-5)
@lru_cache(maxsize=None)
def certified_group(nu_text, ratio_text, count=ROOTS, digits=100):
nu, ratio = QQ(nu_text), QQ(ratio_text)
bits = numberdb.bits(digits,losing=GUARD_BITS)
B = RealBallField(bits)
candidates = approximate_roots(nu,ratio,count)
balls, evidence = [], []
radius = QQ(1)/10**(digits+10)
previous_upper = QQ(0)
for k, guess in enumerate(candidates,1):
low_text = high_precision_root(nu,ratio,guess,digits+30)
high_text = high_precision_root(nu,ratio,guess,digits+60)
low_root, high_root = decimal_rational(low_text), decimal_rational(high_text)
assert abs(low_root-high_root) < max(QQ(1),abs(high_root))/10**(digits+20)
scale = 10**(digits+12)
centre = (high_root*scale).floor()/QQ(scale)
lower, upper = centre-radius, centre+radius
assert previous_upper < lower
s1 = strict_sign(cross_product(nu,ratio,lower,bits))
s2 = strict_sign(cross_product(nu,ratio,upper,bits))
assert s1*s2 == -1
ball = B(centre,radius)
assert ball.is_finite() and ball.accuracy() >= numberdb.bits(digits+8,losing=0)
balls.append(ball)
evidence.append({'k':k,'lower':str(lower),'upper':str(upper),
'endpoint_signs':[s1,s2],
'mpmath_low':low_text,'mpmath_high':high_text,
'accuracy_bits':int(ball.accuracy())})
previous_upper = upper
nodes, certificate = radial_node_count(nu,ratio,previous_upper,bits)
if nodes != count:
raise ArithmeticError('Sturm count %s disagrees with %s proposed roots' % (nodes,count))
return balls, {'nu':nu_text,'lambda':ratio_text,'bits':bits,
'root_count':nodes,'sturm':certificate,'roots':evidence}
class BesselCrossProductZeros(numberdb.Generator):
table = 'T487'
parameters = ('nu','lambda','k')
type = 'R'
digits = 100
rigour = 'proven'
def enumerate(self):
for nu in range(MAX_ORDER+1):
for ratio in RATIOS:
for k in range(1,ROOTS+1):
yield {'nu':str(nu),'lambda':ratio,'k':str(k)}
def value(self,params,digits):
k = int(params['k'])
balls, _ = certified_group(str(params['nu']),str(params['lambda']),max(ROOTS,k),digits)
return balls[k-1]
if __name__ == '__main__':
if os.environ.get('NUMBERDB_KEY_FROM_STDIN') == '1':
numberdb.configure(api_key=sys.stdin.read().strip())
generator = BesselCrossProductZeros()
if '--compute' in sys.argv:
from numberdb._write import to_text
for nu in range(MAX_ORDER+1):
for ratio in RATIOS:
values, proof = certified_group(str(nu),ratio)
proof['written'] = [to_text(v,digits=100) for v in values]
print(json.dumps(proof),flush=True)
elif '--fill' in sys.argv or os.environ.get('NUMBERDB_PUBLISH') == '1':
print(generator.publish(message='Fill the draft with 480 certified Bessel cross-product zeros.'))
else:
report = generator.verify(sample=None)
print(report)
sys.exit(0 if report.ok else 1)