back to table · edit · history · where entries came from · files · download
7183 bytes, as of the version from 2026-10-28 18:59 (current). Recorded here, not run.
"""Independent series checks for kei_0 positive zeros -- numberdb.org/T587.
sage -python checks.py # compute local generator, check series
sage -python checks.py --stored # check every decimal read from the site
For a private draft set NUMBERDB_API_KEY, or NUMBERDB_KEY_FROM_STDIN=1 and pipe
the key. No write occurs. The real series code does not call any Bessel or
Kelvin routine. DLMF 10.65.2 and 10.65.5 are summed with rigorous geometric
bounds on both tails, using H_k <= H_{N+1} + k-N-1 in the harmonic tail.
Every tested ball must be finite and have an absolute radius below 10^-200.
The 100-digit intervals are checked at their written endpoints, not midpoint.
"""
import json
import os
import sys
from decimal import Decimal
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
SERIES_BITS = 1024
TAIL_TARGET = QQ(1)/10**240
RADIUS_LIMIT = QQ(1)/10**200
AS_ROOTS = ('3.91467', '8.34422', '12.78256', '17.22314', '21.66464')
def rational(text):
"""Parse a decimal exactly, without Sage decimal-string coercion."""
d = Decimal(text)
numerator, denominator = d.as_integer_ratio()
return QQ(numerator)/QQ(denominator)
def decimal_interval(text):
text = str(text).rstrip('?')
center = rational(text)
exponent = Decimal(text).as_tuple().exponent
ulp = QQ(10)**exponent
return center-ulp, center+ulp, ulp
def sign(v):
assert v.is_finite() and QQ(v.rad()) < RADIUS_LIMIT
if v > 0:
return 1
if v < 0:
return -1
raise ArithmeticError('Independent series sign uncertain')
def series_parts(t, bits=SERIES_BITS):
R = RealBallField(bits)
t = R(t)
assert t >= 0
term = R(1)
harmonic = R(0)
real, imag, hreal, himag = R(1), R(0), R(0), R(0)
for k in range(1, 3001):
term *= t/(k*k)
harmonic += R(1)/k
residue = k % 4
if residue == 0:
real += term
hreal += harmonic*term
elif residue == 1:
imag += term
himag += harmonic*term
elif residue == 2:
real -= term
hreal -= harmonic*term
else:
imag -= term
himag -= harmonic*term
q = t/((k+2)*(k+2))
if q < QQ(1)/2:
nxt = term*t/((k+1)*(k+1))
tail_i = nxt/(1-q)
tail_h = nxt*((harmonic + R(1)/(k+1))/(1-q) + q/(1-q)**2)
if tail_i + tail_h < TAIL_TARGET:
# RealBall(midpoint, radius) also preserves the original radius.
return (R(real, tail_i.upper()), R(imag, tail_i.upper()),
R(hreal, tail_h.upper()), R(himag, tail_h.upper()), k,
tail_i+tail_h)
raise ArithmeticError('Series tail did not converge')
def kei_series(x, bits=SERIES_BITS):
R = RealBallField(bits)
assert x > 0
ber, bei, hber, hbei, terms, tail = series_parts(R(x)**2/4, bits)
result = -(R(x/2).log()+R.euler_constant())*bei - R.pi()*ber/4 + hbei
assert result.is_finite() and QQ(result.rad()) < RADIUS_LIMIT
return result, terms, tail
def controls():
R = RealBallField(SERIES_BITS)
C = ComplexBallField(SERIES_BITS)
ber, bei, hber, hbei, _, _ = series_parts(0)
assert ber == 1 and bei == 0 and hber == 0 and hbei == 0
# A nonzero closed form tests the special-function engine used by the builder.
khalf = C(1).bessel_K(QQ(1)/2).real()
known = (R.pi()/2).sqrt()*(-R(1)).exp()
assert khalf.is_finite() and known.is_finite()
assert QQ(khalf.rad()) < RADIUS_LIMIT and QQ(known.rad()) < RADIUS_LIMIT
assert khalf.overlaps(known) and khalf > QQ(2)/5 and khalf < QQ(1)/2
x = QQ(1)/2**128
v, _, _ = kei_series(x)
# The defining limit is -pi/4; the error at this exact x is O(x^2 |log x|).
assert abs(v+R.pi()/4) < QQ(1)/10**73
assert sign(v) == -1
return {'series_at_origin': 'I=1, harmonic sum=0',
'nonzero_control': 'K_(1/2)(1)=sqrt(pi/2)/e, finite tight overlap',
'origin_limit': 'kei_0(2^-128) differs from -pi/4 by less than 1e-73'}
def check_rows(rows):
control = controls()
assert [int(r['n']) for r in rows] == list(range(1, 21))
evidence = []
previous = QQ(0)
for row in rows:
n = int(row['n'])
left, right, ulp = decimal_interval(row['number'])
assert left > previous
previous = right
fl, kl, tl = kei_series(left)
fr, kr, tr = kei_series(right)
sl, sr = sign(fl), sign(fr)
assert sl == (-1)**n and sr == (-1)**(n+1)
# Move the entire interval ten last-place units to the right.
nl, _, _ = kei_series(left+10*ulp)
nr, _, _ = kei_series(right+10*ulp)
assert sign(nl) == sign(nr) == sr
as_check = None
if n <= 5:
a, b, _ = decimal_interval(AS_ROOTS[n-1])
assert a < left < right < b
as_check = AS_ROOTS[n-1]
# Check the differential identity and its phase sign at each row.
R = RealBallField(SERIES_BITS)
C = ComplexBallField(SERIES_BITS)
center = (left+right)/2
a = C(1,1)/R(2).sqrt()
z = C(center)*a
g = z.bessel_K(0)
gp = -a*z.bessel_K(1)
gpp = a*a*(z.bessel_K(0)+z.bessel_K(2))/2
residual = gpp+gp/center-C(0,1)*g
for component in (residual.real(),residual.imag()):
assert component.is_finite() and component.contains_zero()
assert QQ(component.rad()) < RADIUS_LIMIT
w = g.real()*gp.imag()-g.imag()*gp.real()
assert w.is_finite() and w < 0
evidence.append({'n':n,'number':str(row['number']),
'ode_identity_contains_zero':True,
'phase_derivative_negative':True,
'series_endpoint_signs':[sl,sr],
'series_terms':max(kl,kr),
'series_max_radius':str(max(fl.rad(),fr.rad())),
'series_tail_sum_bound':str(max(tl.upper(),tr.upper())),
'displaced_interval_same_sign':True,
'abramowitz_stegun_5dp':as_check})
return {'controls':control,'independent_series_bits':SERIES_BITS,
'series_radius_limit':'1e-200','rows':evidence,
'all_20_written_intervals_certified':True}
if __name__ == '__main__':
if os.environ.get('NUMBERDB_KEY_FROM_STDIN') == '1':
numberdb.configure(api_key=sys.stdin.read().strip())
if '--stored' in sys.argv:
d=numberdb.table('T587')
block=d['Numbers']
if isinstance(block,list):
rows=[{'n':int(r['params']['n']),'number':r['number']} for r in block]
else:
rows=[{'n':int(k),'number':v['number']} for k,v in block.items()]
else:
import generate
from numberdb._write import to_text
g=generate.KelvinKeiZeros()
rows=[{'n':p['n'],'number':to_text(g.value(p,100),digits=100)}
for p in g.enumerate()]
print(json.dumps(check_rows(rows),indent=2))