back to table · edit · history · where entries came from · files · download
5056 bytes, as of the version from 2026-10-28 20:22 (current). Recorded here, not run.
"""Positive zeros of the Kelvin function ber_0 -- numberdb.org/T589.
Run with SageMath (install numberdb once with sage -pip install numberdb):
sage -python generate.py # verify all stored entries, read only
sage -python checks.py # independent offline certificates
KelvinBerZeros().value({'n': 21}, 100) extends the same indexed family.
No command in this file creates, offers, reviews, or publishes a table.
All precision claims come from Arb endpoint signs and the complete root-count
certificate, not from the mpmath candidate finder. See certificate.md and
argument-selection.md. 800 working bits retain ample precision for the 100
written digits through n=20; the narrow brackets have radius 10^(-digits-12).
"""
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
from fractions import Fraction
from functools import lru_cache
import mpmath as mp
WORKING_BITS = 800
DEFAULT_COUNT = 20
def rational_decimal(text):
f = Fraction(text)
return QQ(f.numerator) / QQ(f.denominator)
def sign(v):
if not v.is_finite():
raise ArithmeticError('nonfinite function enclosure')
if v > 0:
return 1
if v < 0:
return -1
raise ArithmeticError('function sign is not certified')
def fields(digits):
bits = max(WORKING_BITS, numberdb.bits(digits, losing=400))
return RealBallField(bits), ComplexBallField(bits)
def solution(x, R, C):
c = C(1, 1) / R(2).sqrt()
z = c * C(x)
return z.bessel_I(0), c * z.bessel_I(1)
@lru_cache(maxsize=None)
def counted_brackets(count=DEFAULT_COUNT, bits=WORKING_BITS):
"""Exhaustively isolate the first count roots, with no skipped cells.
F'' = i F - F'/x. Thus (x W)' = x |F|^2 for
W = Re(F) Im(F') - Im(F) Re(F'), and W>0 on x>0.
The continuous phase strictly increases. On [a,a+h], a>=1,
||(F,F')||_infinity <= N exp((1+1/a)h), N=|F(a)|+|F'(a)|.
Consequently |F(t)/F(a)-1| <= h N exp((1+1/a)h)/|F(a)|.
Proving this less than 1 restricts the phase change to less than pi/2,
so each cell contains at most one zero of Re(F). Endpoint signs then
decide its exact root count. On [0,1], |F-1|<=exp(1/4)-1<1.
"""
if count < 1:
raise ValueError('root index must be positive')
R, C = RealBallField(bits), ComplexBallField(bits)
assert (R(QQ(1)/4).exp()-1) < 1
assert C(0).bessel_I(0) == 1 # nonzero known-answer control
a, h = QQ(1), QQ(1)/8
fa, fpa = solution(a, R, C)
sa = sign(fa.real())
assert sa == 1
brackets, certificates = [], []
while len(brackets) < count:
b = a+h
bound = R(h)*(R((1+1/a)*h)).exp()*(abs(fa)+abs(fpa))/abs(fa)
assert bound.is_finite() and bound < 1
fb, fpb = solution(b, R, C)
sb = sign(fb.real())
certificates.append({'a': str(a), 'b': str(b),
'relative_disk_bound': str(bound),
'left_sign': sa, 'right_sign': sb})
if sa != sb:
brackets.append((a,b))
a, fa, fpa, sa = b, fb, fpb, sb
return tuple(brackets), tuple(certificates)
@lru_cache(maxsize=None)
def certified_root(n, digits=100):
R, C = fields(digits)
brackets, _ = counted_brackets(max(DEFAULT_COUNT,n), R.precision())
a, b = brackets[n-1]
with mp.workdps(digits+50):
left = mp.mpf(int(a.numerator()))/int(a.denominator())
right = mp.mpf(int(b.numerator()))/int(b.denominator())
candidate = mp.findroot(lambda x: mp.ber(0,x), (left,right))
centre = rational_decimal(mp.nstr(candidate, digits+35))
radius = QQ(1)/10**(digits+12)
low, high = centre-radius, centre+radius
assert a < low < high < b
sl = sign(solution(low,R,C)[0].real())
sh = sign(solution(high,R,C)[0].real())
assert sl == (-1)**(n-1) and sh == -sl
enclosure = RealIntervalField(R.precision())(low,high)
assert enclosure.absolute_diameter() > 0
return enclosure
class KelvinBerZeros(numberdb.Generator):
table = 'T589'
parameters = ('n',)
type = 'R'
digits = 100
rigour = 'proven'
files = ('generate.py', 'checks.py', 'certificate.md', 'argument-selection.md',
'skill-evidence.json')
def enumerate(self, count=DEFAULT_COUNT):
for n in range(1,count+1):
yield {'n': n}
def value(self, params, digits):
n = int(params['n'])
if n < 1:
raise ValueError('root index must be positive')
return certified_root(n,digits)
if __name__ == '__main__':
import os, sys
if os.environ.get('NUMBERDB_KEY_FROM_STDIN') == '1':
key = sys.stdin.read().strip()
if '=' in key:
key = key.split('=',1)[1].strip().strip('\"\'')
os.environ['NUMBERDB_API_KEY'] = key
report = KelvinBerZeros().verify(sample=None)
print(report)
sys.exit(0 if report.ok else 1)