back to table · edit · history · where entries came from · files · download
5156 bytes, as of the version from 2026-10-28 18:40. Recorded here, not run.
"""Positive zeros of the Kelvin function kei_0 -- numberdb.org/T587.
Run with SageMath or passagemath and numberdb installed:
sage -pip install numberdb # once
sage -python generate.py # verify every stored entry
sage -python generate.py --fill # write entries, key in environment
sage -python checks.py # independent series certificates
The table has every root index n=1,...,20. See argument-selection.md for the
extent and root-counting.md for the proof that no positive roots were skipped.
The first strictly positive root has index 1. G(x)=K_0(x*(1+i)/sqrt(2)), with
principal K_0, and the stored number is x, not x/sqrt(2).
A 1/4 mesh is used only to certify the count, never as table arguments.
The image of each cell is certified to lie in a disk excluding zero using
|G'(x)| <= K_1(x/sqrt(2)). Strict phase monotonicity then makes the sign count
complete. Newton iterates merely locate candidates. Exact rational endpoints
with proven opposite signs certify each returned root enclosure.
Working precision is bits(digits, losing=420); the counting scan uses 384 bits.
Failure of any bound raises rather than guessing. The measured 100-digit run
is documented in numerical-evidence.json and the independent checks use 1024
bits with explicit series tails. All source dependencies are attached.
"""
import os
import sys
from functools import lru_cache
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
COUNT_BITS = 384
WORKING_GUARD = 420
DEFAULT_LIMIT = 20
def finite(z):
if hasattr(z, 'imag'):
assert z.real().is_finite() and z.imag().is_finite()
else:
assert z.is_finite()
return z
def sign(v):
finite(v)
if v > 0:
return 1
if v < 0:
return -1
raise ArithmeticError('Function sign is not certified')
def field_pair(bits):
R = RealBallField(bits)
C = ComplexBallField(bits)
return R, C, C(1, 1) / R(2).sqrt()
def kei(x, bits):
R, C, a = field_pair(bits)
return finite((C(x) * a).bessel_K(0)).imag()
@lru_cache(maxsize=None)
def counted_brackets(limit=DEFAULT_LIMIT):
"""Certify every cell from 1/4 through the requested root."""
assert isinstance(limit, int) and limit >= 1
R, C, a = field_pair(COUNT_BITS)
# Analytic origin bound in root-counting.md, checked exactly here.
assert QQ(5)/63 + QQ(64)/3969 < QQ(3)/4
l = QQ(1)/4
step = QQ(1)/4
left_sign = sign((C(l)*a).bessel_K(0).imag())
assert left_sign == -1
brackets = []
cells = []
while len(brackets) < limit:
r = l + step
m = (l + r)/2
g = finite((C(m)*a).bessel_K(0))
bound = finite(C(R(l)/R(2).sqrt()).bessel_K(1)).real()
radius = R(step/2)*bound
assert radius > 0 and radius < abs(g), ('cell disk failed', l, r)
right_sign = sign((C(r)*a).bessel_K(0).imag())
cells.append((l, r, radius/abs(g), left_sign, right_sign))
if right_sign != left_sign:
brackets.append((l, r))
l, left_sign = r, right_sign
return tuple(brackets), tuple(cells)
def root_enclosure(n, digits=100):
n = int(n)
assert n >= 1
brackets, cells = counted_brackets(max(n, DEFAULT_LIMIT))
lo, hi = brackets[n-1]
bits = numberdb.bits(digits, losing=WORKING_GUARD)
R, C, a = field_pair(bits)
m = (lo+hi)/2
target = QQ(1)/10**(digits+12)
for iteration in range(20):
z = C(m)*a
f = finite(z.bessel_K(0)).imag()
derivative = finite(-a*z.bessel_K(1)).imag()
assert not derivative.contains_zero()
correction = finite(f/derivative)
m -= correction.mid().exact_rational()
assert lo < m < hi
if abs(correction) < target:
break
else:
raise ArithmeticError('Newton candidate did not stabilize')
eps = QQ(1)/10**(digits+8)
l, r = m-eps, m+eps
assert lo < l < r < hi
assert sign(kei(l, bits))*sign(kei(r, bits)) == -1
value = finite(R(m, eps))
assert value.rad() > 0 and QQ(value.rad()) < 2*eps
return value
class KelvinKeiZeros(numberdb.Generator):
table = 'T587'
parameters = ('n',)
type = 'R'
digits = 100
rigour = 'proven'
files = ('generate.py', 'checks.py', 'argument-selection.md',
'root-counting.md', 'skill-evidence.json', 'numerical-evidence.json')
def enumerate(self, limit=DEFAULT_LIMIT):
for n in range(1, int(limit)+1):
yield {'n': n}
def value(self, params, digits):
return root_enclosure(int(params['n']), digits)
if __name__ == '__main__':
if os.environ.get('NUMBERDB_KEY_FROM_STDIN') == '1':
numberdb.configure(api_key=sys.stdin.read().strip())
gen = KelvinKeiZeros()
if '--fill' in sys.argv or os.environ.get('NUMBERDB_PUBLISH') == '1':
print(gen.publish(message='Fill twenty certified positive kei_0 roots'))
else:
report = gen.verify(sample=None)
print(report)
sys.exit(0 if report.ok else 1)