back to table · edit · history · where entries came from · files · download
5275 bytes, as of the version from 2026-10-28 19:40 (current). Recorded here, not run.
"""Positive zeros of the Kelvin function ker_0 -- numberdb.org/T588.
Run with SageMath (or a suitable passagemath installation):
sage -pip install numberdb mpmath
sage -python generate.py # verify every stored entry
sage -python generate.py --fill # write entries with NUMBERDB_API_KEY
The --fill option fills the existing table; it never offers or publishes it.
For an unpublished draft, verification also needs its author's API key.
The principal-branch convention is ker_0(x)=Re K_0(x exp(pi i/4)). The
zero index starts at one and the default range is 1,...,20. See
argument-selection.md and root-counting.md for range and certificates.
384-bit ball computations certify a complete scan by an origin bound,
phase monotonicity and disks excluding zero on consecutive cells. mpmath
only proposes a narrow bracket. Arb certifies its endpoint signs with a
420-bit guard beyond the requested digits (753 bits at 100 digits).
The returned interval contains the rational bracket, radius 10^(-digits-8).
The independent checks.py evaluates explicit real series with tail bounds.
"""
import os
import sys
from functools import lru_cache
from fractions import Fraction
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 sage.rings.real_mpfi import RealIntervalField
import mpmath as mp
COUNT_BITS = 384
WORKING_GUARD = 420
DEFAULT_LIMIT = 20
def rational_decimal(text):
value = Fraction(str(text))
return QQ(value.numerator) / value.denominator
def sign(value):
assert value.is_finite(), 'nonfinite enclosure'
if value > 0:
return 1
if value < 0:
return -1
raise ArithmeticError('sign not certified')
def fields(bits):
R = RealBallField(bits)
C = ComplexBallField(bits)
return R, C, C(1, 1) / R(2).sqrt()
@lru_cache(maxsize=None)
def counted_cells(limit=DEFAULT_LIMIT):
"""Return one certified cell for each root and complete-count evidence."""
R, C, s = fields(COUNT_BITS)
t = QQ(1) / 64
A = t / (1-t)
B = t / (1-t)**2
origin_lower = (R(8).log()-R.euler_constant())*(1-A)-R.pi()*A/4-B
assert origin_lower > 0
cells = []
roots = []
a = QQ(1)/4
ga = (s*a).bessel_K(0)
sa = sign(ga.real())
assert sa == 1
while len(roots) < limit:
b = a+QQ(1)/4
for attempt in range(30):
m = (a+b)/2
gm = (s*m).bessel_K(0)
derivative_bound = (C(a)/R(2).sqrt()).bessel_K(1).real()
radius = ((b-a)/2)*derivative_bound
if (gm.real().is_finite() and gm.imag().is_finite()) and radius < abs(gm):
break
b = m
else:
raise ArithmeticError('could not certify scan cell')
gb = (s*b).bessel_K(0)
sb = sign(gb.real())
ratio = radius/abs(gm)
assert ratio < 1
cells.append({'a':str(a),'b':str(b),'left_sign':sa,'right_sign':sb,
'disk_ratio_upper':str(ratio.upper())})
if sa != sb:
roots.append((a,b))
a, ga, sa = b, gb, sb
assert len(cells) < 10000*limit
return tuple(roots), {'origin_lower':str(origin_lower), 'cells':cells,
'end':str(a), 'count':len(roots)}
@lru_cache(maxsize=None)
def enclosed_root(n, digits=100):
n = int(n)
assert n >= 1
brackets, _ = counted_cells(max(DEFAULT_LIMIT, n))
a, b = brackets[n-1]
with mp.workdps(digits+35):
ma = mp.mpf(int(a.numerator()))/int(a.denominator())
mb = mp.mpf(int(b.numerator()))/int(b.denominator())
proposal = mp.findroot(lambda x: mp.ker(0,x), (ma,mb))
center = rational_decimal(mp.nstr(proposal,digits+25))
epsilon = QQ(1)/10**(digits+8)
lo, hi = center-epsilon, center+epsilon
assert a < lo < hi < b
bits = numberdb.bits(digits,losing=WORKING_GUARD)
R, C, s = fields(bits)
left = (s*lo).bessel_K(0).real()
right = (s*hi).bessel_K(0).real()
assert sign(left)*sign(right) == -1
return RealIntervalField(bits)(lo,hi)
class KelvinKerZeros(numberdb.Generator):
table = 'T588'
parameters = ('n',)
type = 'R'
digits = 100
rigour = 'proven'
files = ('generate.py','checks.py','argument-selection.md','root-counting.md')
def enumerate(self, limit=DEFAULT_LIMIT):
for n in range(1,limit+1):
yield {'n':n}
def value(self, params, digits):
return enclosed_root(int(params['n']),digits)
if __name__ == '__main__':
if os.environ.get('NUMBERDB_KEY_FROM_STDIN') == '1':
raw = sys.stdin.read().strip()
if '=' in raw and raw.split('=',1)[0].strip().isidentifier():
raw = raw.split('=',1)[1].strip().strip('\"\'')
numberdb.configure(api_key=raw)
generator = KelvinKerZeros()
if '--fill' in sys.argv or os.environ.get('NUMBERDB_PUBLISH') == '1':
print(generator.publish(message='Fill twenty certified positive ker_0 zeros',
assisted_by=os.environ.get('NUMBERDB_ASSISTED_BY','')))
else:
report = generator.verify(sample=None)
print(report)
sys.exit(0 if report.ok else 1)