back to table · edit · history · where entries came from · files · download
9746 bytes, as of the version from 2026-10-25 15:06. Recorded here, not run.
"""Independent checks for Kelvin bei -- numberdb.org/T491.
sage -python checks.py
sage -python checks.py --stored # NUMBERDB_API_KEY for the draft
No writes. Requires generate.py and check.py beside this file, and mpmath.
An exact-rational series with a geometric tail bound is independent of Arb's
Bessel evaluator. Bessel's integral is checked by numerical quadrature at
125 decimal digits, with nonzero and zero controls first. The quadrature and
mpmath comparisons are independent numerical checks, not error certificates.
The certified series and the generator balls establish interval containment.
DLMF 10.65.1: a_k = (x/2)^n (x^2/4)^k / (k! (n+k)!).
For k >= N, a_(k+1)/a_k <= r=(x^2/4)/((N+1)(n+N+1)).
Thus the omitted absolute tail is <= a_N/(1-r), since |sin|,|cos| <= 1.
Trigonometric factors at multiples of pi/4 use an exact sign cycle and,
only for odd n, a 600-bit enclosure of 1/sqrt(2).
"""
import json
import math
import os
import sys
from fractions import Fraction
from pathlib import Path
import numberdb.sage as numberdb
from sage.rings.complex_arb import ComplexBallField
from sage.rings.rational_field import QQ
from sage.rings.real_arb import RealBallField
from sage.env import SAGE_VERSION
import mpmath as mp
from numberdb._write import to_text
from generate import KelvinBei
from check import exactness, measure, names_its_rings
R = RealBallField(600)
C = ComplexBallField(600)
NTERMS = 200
def rational_decimal(text):
q = Fraction(str(text))
return QQ(q.numerator) / QQ(q.denominator)
def written_bounds(text):
"""Exact endpoints of NumberDB's decimal +/- one last-place unit."""
text = str(text)
mantissa, _, exponent = text.lower().partition('e')
decimal_places = len(mantissa.split('.')[1]) if '.' in mantissa else 0
ulp = QQ(10) ** ((int(exponent) if exponent else 0) - decimal_places)
center = rational_decimal(text)
return center - ulp, center + ulp
def as_ball(text):
lo, hi = written_bounds(text)
return R((lo + hi) / 2).add_error(R((hi - lo) / 2))
def encloses(text, value):
assert value.is_finite()
lo, hi = written_bounds(text)
assert R(lo) < value and value < R(hi), (text, str(value))
def series_pair(n, x):
"""Certified ber, bei enclosures, without calling any special function."""
x = QQ(x)
a = (x / QQ(2))**n / QQ(math.factorial(n))
t = x*x / QQ(4)
real = QQ(0)
imag = QQ(0)
sin_sign = (0, 1, 1, 1, 0, -1, -1, -1)
cos_sign = (1, 1, 0, -1, -1, -1, 0, 1)
for k in range(NTERMS):
phase = (3*n + 2*k) % 8
real += cos_sign[phase] * a
imag += sin_sign[phase] * a
a *= t / QQ((k+1)*(n+k+1))
ratio = t / QQ((NTERMS+1)*(n+NTERMS+1))
assert 0 <= ratio < 1
tail = a / (1-ratio)
factor = 1 / R(2).sqrt() if n % 2 else R(1)
return ((R(real)*factor).add_error(R(tail)),
(R(imag)*factor).add_error(R(tail)), tail)
def bessel_integral(n, z):
"""DLMF 10.9.2, a quadrature with elementary integrand only."""
return mp.quad(lambda theta: mp.cos(z*mp.sin(theta)-n*theta),
[0, mp.pi/2, mp.pi]) / mp.pi
def leaf_value(leaf):
return str(leaf['number']) if isinstance(leaf, dict) else str(leaf)
def main(stored=False):
if os.environ.get('NUMBERDB_KEY_FROM_STDIN') == '1':
os.environ['NUMBERDB_API_KEY'] = sys.stdin.read().strip()
g = KelvinBei()
params = list(g.enumerate())
assert len(params) == 360
assert not names_its_rings(str(Path(__file__).with_name('generate.py')))
values = {(p['nu'], p['x']):g.value(p, g.digits) for p in params}
assert not exactness(values)
texts = {key:to_text(v, digits=100) for key,v in values.items()}
ber_texts = None
if stored:
doc = numberdb.table(g.table)
assert list(doc['Parameters']) == ['nu','x']
assert doc['Data properties']['digits'] == '100'
assert doc['Data properties']['rigour'] == 'proven'
texts = {(n,x):leaf_value(v) for n,group in doc['Numbers'].items()
for x,v in group.items()}
assert set(texts) == set(values)
ber = numberdb.table('T490')
ber_texts = {(n,x):leaf_value(v) for n,group in ber['Numbers'].items()
for x,v in group.items()}
assert set(ber_texts) == set(values)
# A nonzero control detects an accidentally zero-valued measurement.
s0 = series_pair(0, QQ(0))
assert s0[0] == 1 and s0[1] == 0 and s0[2] == 0
assert series_pair(1, QQ(0))[1] == 0
mp.mp.dps = 125
assert abs(bessel_integral(0, mp.mpf(0)) - 1) < mp.mpf('1e-120')
for n in (1,2):
assert abs(bessel_integral(n, mp.mpf(0))) < mp.mpf('1e-120')
print('Controls passed: rational series (ber, bei)=(1,0), integrals J_0(0)=1 and J_1(0)=J_2(0)=0.',flush=True)
series = {}
max_radius = R(0)
max_relative_radius = R(0)
max_tail = QQ(0)
worst_quad_relative = mp.mpf(0)
worst_mp_relative = mp.mpf(0)
for i,p in enumerate(params,1):
key = (p['nu'],p['x'])
n,x = int(p['nu']),QQ(p['x'])
v, text = values[key], texts[key]
assert v.is_finite() and not v.contains_zero()
relative = R(v.rad()) / abs(R(v))
assert relative < R(QQ(1)/10**110)
max_radius = max(max_radius, R(v.rad()))
max_relative_radius = max(max_relative_radius, relative)
encloses(text, R(v))
real, imaginary, tail = series_pair(n,x)
assert real.is_finite() and imaginary.is_finite()
assert R(imaginary.rad()) < abs(imaginary)*R(QQ(1)/10**150)
assert tail < QQ(1)/10**150
max_tail = max(max_tail,tail)
assert R(v).overlaps(imaginary)
encloses(text, imaginary)
series[key] = (real,imaginary)
if ber_texts is not None:
encloses(ber_texts[key], real)
# Independent modified Bessel representation, DLMF 10.61.1.
rotated = (C(0,1)**n * (C(x)*C(1,1)/C(2).sqrt()).bessel_I(C(n))).imag()
assert rotated.is_finite() and R(rotated.rad()) < abs(rotated)*R(QQ(1)/10**150)
encloses(text,rotated)
# Integer-order reflection identities, checked on every row.
z = C(x)*C(-1,1)/C(2).sqrt()
reflected_order = ((-1)**n)*z.bessel_J(C(-n)).imag()
reflected_argument = ((-1)**n)*(-z).bessel_J(C(n)).imag()
encloses(text,reflected_order)
encloses(text,reflected_argument)
with mp.workdps(160):
mx = mp.mpf(int(x.numerator())) / int(x.denominator())
mv = mp.bei(n,mx)
point = rational_decimal(mp.nstr(mv,155))
lo,hi = written_bounds(text)
assert lo < point < hi
center = mp.mpf(str(v.mid()))
mp_error = abs(mv-center)/abs(mv)
assert mp_error < mp.mpf('1e-110')
worst_mp_relative = max(worst_mp_relative,mp_error)
mx = mp.mpf(int(x.numerator())) / int(x.denominator())
zq = mx * mp.mpc(-1,1) / mp.sqrt(2)
integral = mp.im(bessel_integral(n,zq))
assert mp.isfinite(integral)
quad_error = abs(integral-mv)/abs(mv)
assert quad_error < mp.mpf('1e-110'), (key,quad_error)
worst_quad_relative = max(worst_quad_relative,quad_error)
lo,hi = written_bounds(text)
assert lo < rational_decimal(mp.nstr(integral,120)) < hi
if i % 60 == 0:
print('checked %d/360: series, I-rotation, reflections, mpmath, Bessel integral' % i,flush=True)
# Recurrence with all three orders present, on written intervals.
arguments = [p['x'] for p in params if p['nu']=='0']
max_recurrence_scaled_radius = R(0)
for xs in arguments:
x=QQ(xs)
b0=as_ball(texts[('0',xs)])
b1=as_ball(texts[('1',xs)])
b2=as_ball(texts[('2',xs)])
ber1=(as_ball(ber_texts[('1',xs)]) if ber_texts is not None
else series[('1',xs)][0])
residual=b0+b2+R(2).sqrt()/R(x)*(b1+ber1)
scale=abs(b0)+abs(b2)+R(2).sqrt()/R(x)*(abs(b1)+abs(ber1))
assert residual.is_finite() and residual.contains_zero()
assert R(residual.rad()) < scale*R(QQ(1)/10**95)
max_recurrence_scaled_radius=max(max_recurrence_scaled_radius,R(residual.rad())/scale)
# Both Programs snippets, as written in the document, are executable.
document_path=Path(__file__).with_name('document.json')
if document_path.exists():
doc=json.loads(document_path.read_text())
for name,program in doc['Programs'].items():
ns={}
exec(program['code'],ns)
result=ns['value']
if name=='program-sage':
encloses(texts[('0','1')],R(result))
else:
lo,hi=written_bounds(texts[('0','1')])
assert lo < rational_decimal(mp.nstr(result,100)) < hi
print(json.dumps({'mode':'stored' if stored else 'local',
'table':'T491','entries':len(texts), 'size':measure(texts),
'precision_bits':numberdb.bits(100,losing=64),
'maximum_generator_radius':str(max_radius),
'maximum_generator_relative_radius':str(max_relative_radius),
'maximum_series_tail_bound':str(R(max_tail)),
'maximum_mpmath_relative_error':mp.nstr(worst_mp_relative,10),
'maximum_quadrature_relative_error':mp.nstr(worst_quad_relative,10),
'recurrences_checked':len(arguments),
'maximum_recurrence_scaled_radius':str(max_recurrence_scaled_radius),
'Sage':SAGE_VERSION, 'mpmath':mp.__version__,
'numberdb':numberdb.__version__ if hasattr(numberdb,'__version__') else 'see package metadata',
'passed':True},indent=2),flush=True)
return texts
if __name__=='__main__':
main(stored='--stored' in sys.argv)