back to table · edit · history · where entries came from · files · download
10179 bytes, as of the version from 2026-10-25 16:35. Recorded here, not run.
"""Independent numerical checks for numberdb.org/T494.
sage -python checks.py
sage -python checks.py --stored [--key-stdin]
Requires generate.py and the NumberDB table-build check.py alongside this file.
The offline run prints the written values and independent integral enclosures
as EVIDENCE_JSON, and a compact SUMMARY_JSON. --stored repeats all comparisons
on the decimal intervals read from T494, using stored I_0,I_1 from T196 where
its argument range overlaps. No writes are performed.
The principal independent test is certified 480-bit quadrature of DLMF 11.5.6.
At these orders 2*nu is a nonnegative integer, so the integrand is entire.
A separate 600-bit calculation integrates the exact rational series of I_0
with a geometric tail bound. This checks the cross-table integral identity.
mpmath at 160 digits provides another implementation, not an error certificate.
"""
import os
import sys
import json
from decimal import Decimal
import numberdb.sage as numberdb
from sage.rings.rational_field import QQ
from sage.rings.complex_arb import ComplexBallField
from sage.rings.real_arb import RealBallField
from sage.rings.real_mpfi import RealIntervalField
import sage.version
import mpmath as mp
from numberdb._write import to_text
from numberdb._generate import _stored_entries
from numberdb._compare import digits_of
from generate import ModifiedStruve, struve_l
from check import exactness, measure, names_its_rings
C = ComplexBallField(600)
R = RealBallField(600)
CI = ComplexBallField(480)
RI = RealIntervalField(600)
mp.mp.dps = 160
def as_rational_decimal(text):
a, b = Decimal(text).as_integer_ratio()
return QQ(a) / QQ(b)
def bounds(text):
d = Decimal(text)
m = as_rational_decimal(text)
r = QQ(10) ** int(d.as_tuple().exponent)
return m - r, m + r
def written_ball(text):
lo, hi = bounds(text)
return R((lo + hi) / 2, (hi - lo) / 2)
def sane(v, relative_digits=110):
assert v.is_finite(), str(v)
assert not v.contains_zero(), str(v)
assert v.rad_as_ball() < abs(v) * R(10) ** (-relative_digits), str(v)
def real(v):
assert v.real().is_finite() and v.imag().is_finite()
assert v.imag().contains_zero(), str(v)
return v.real()
def inside(text, value):
lo, hi = bounds(text)
interval = RI(value)
return lo <= interval.lower().exact_rational() and interval.upper().exact_rational() <= hi
def integral_l(nu, x):
assert 2 * nu in (0, 1, 2, 3, 4)
z = CI(x)
coefficient = 2 * (z / 2) ** CI(nu) / (CI.pi().sqrt() * CI(nu + QQ(1)/2).gamma())
value = coefficient * CI.integral(
lambda theta, analytic: (z * theta.cos()).sinh() * theta.sin() ** int(2 * nu),
0, CI.pi() / 2)
return real(value)
def half_order(nu, x):
z = R(x)
factor = (2 / (R.pi() * z)).sqrt()
if nu == QQ(1)/2:
return factor * (z.cosh() - 1)
assert nu == QQ(3)/2
return factor * (z.sinh() - (z.cosh() - 1) / z - z / 2)
def integral_i0(x):
# Terms x^(2k+1)/(4^k (k!)^2 (2k+1)), k=0,...,249.
x = QQ(x)
term = x
total = QQ(0)
for k in range(250):
total += term
term *= (x*x/4) * QQ(2*k+1) / ((k+1)**2 * (2*k+3))
# term is the first omitted term. Subsequent ratios are bounded by r.
ratio = (x*x/4) / QQ(251**2)
assert ratio < 1
tail = term / (1-ratio)
return R(total, tail), tail
def residual_ok(v, scale):
assert v.is_finite() and v.contains_zero(), str(v)
assert v.rad_as_ball() < max(R(1), abs(scale)) * R(10)**(-75), str(v)
def controls():
for f, a, b, known in [
(lambda t, analytic: t, 0, 1, QQ(1)/2),
(lambda t, analytic: t.sin(), 0, CI.pi()/2, 1),
]:
value = real(CI.integral(f, a, b))
assert value.contains_exact(known)
sane(value, 130)
zero = integral_l(QQ(0), QQ(0))
assert zero.contains_zero() and zero.rad_as_ball() < R(10)**(-130)
elementary = half_order(QQ(1)/2, QQ(1))
quadrature = integral_l(QQ(1)/2, QQ(1))
sane(quadrature, 130)
assert quadrature.overlaps(elementary)
assert not quadrature.overlaps(elementary + R('0.01'))
zero_i, tail = integral_i0(0)
assert zero_i.is_zero() and tail == 0
independent_i = real(C.integral(lambda t, analytic: t.bessel_I(0), 0, 1))
sane(independent_i, 150)
series_i, _ = integral_i0(1)
assert independent_i.overlaps(series_i)
# Wrong sign and missing pi/2 from the proposal both fail this known test.
l0, l1 = integral_l(QQ(0), QQ(1)), integral_l(QQ(1), QQ(1))
i0, i1 = real(C(1).bessel_I(0)), real(C(1).bessel_I(1))
wrong = i0 * (2/R.pi() + l1) + i1*l0
assert not wrong.overlaps(independent_i)
fixed = R.pi()/2 * (i0*(2/R.pi()+l1) - i1*l0)
residual_ok(fixed-independent_i, independent_i)
print('CONTROLS passed: polynomial and sine integrals, origin, elementary half-order, I0 integral, rejected perturbed value and incorrect proposal identity', flush=True)
def main(stored=False):
controls()
gen = ModifiedStruve()
assert not names_its_rings(os.path.join(os.path.dirname(__file__), 'generate.py'))
parameters = list(gen.enumerate())
assert len(parameters) == 600
keys = [(QQ(p['nu']), QQ(p['x'])) for p in parameters]
assert len(set(keys)) == 600
computed = {key: gen.value(p,100) for key,p in zip(keys,parameters)}
assert not exactness(computed)
text = {key:to_text(value,100) for key,value in computed.items()}
if stored:
doc = numberdb.table(gen.table)
raw = _stored_entries(doc)
text = {tuple(QQ(x) for x in identity.split(',')): value for identity,value in raw.items()}
assert set(text) == set(keys)
values = {key:written_ball(t) for key,t in text.items()}
maxima = {'generator_radius':R(0), 'generator_relative_radius':R(0),
'quadrature_radius':R(0), 'quadrature_relative_radius':R(0),
'i0_tail_bound':R(0)}
counts = {'rows':0,'integral_enclosures':0,'mpmath':0,'half_orders':0,
'recurrences':0,'rotation':0,'negative_order':0,'bessel_integral':0,
'stored_bessel_arguments':0}
evidence = []
for index, key in enumerate(keys):
nu, x = key
v, t = computed[key], text[key]
sane(v)
assert digits_of(t) == 100
assert inside(t,v), key
maxima['generator_radius'] = max(maxima['generator_radius'],R(v.rad()))
maxima['generator_relative_radius'] = max(maxima['generator_relative_radius'],R(v.rad())/abs(v))
integral = integral_l(nu,x)
sane(integral,130)
assert inside(t,integral), ('quadrature',key)
maxima['quadrature_radius'] = max(maxima['quadrature_radius'],R(integral.rad()))
maxima['quadrature_relative_radius'] = max(maxima['quadrature_relative_radius'],R(integral.rad())/abs(integral))
counts['integral_enclosures'] += 1
mn = mp.mpf(int(nu.numerator())) / int(nu.denominator())
mx = mp.mpf(int(x.numerator())) / int(x.denominator())
expected = as_rational_decimal(mp.nstr(mp.struvel(mn,mx),160))
lo,hi = bounds(t)
assert lo < expected < hi, ('mpmath',key)
counts['mpmath'] += 1
if nu.denominator() == 2:
elementary = half_order(nu,x)
sane(elementary,150)
assert inside(t,elementary), ('half-order',key)
counts['half_orders'] += 1
def at(order):
return values.get((order,x), struve_l(order,x,600))
forcing = (R(x)/2)**R(nu)/(R.pi().sqrt()*R(nu+QQ(3)/2).gamma())
lhs = at(nu-1)-at(nu+1)
rhs = (2*nu/x)*values[key]+forcing
residual_ok(lhs-rhs,rhs)
counts['recurrences'] += 1
# DLMF 11.2.1 and 11.2.2, principal branches at i*x.
z, n = C(0,x), C(nu)
h = ((z/2)**(n+1) * (-z*z/4).hypergeometric([1],[QQ(3)/2,n+QQ(3)/2])
/ (C(QQ(3)/2).gamma()*(n+QQ(3)/2).gamma()))
rotated = real(-C(0, 1)*(-C(0, 1)*C.pi()*n/2).exp()*h)
sane(rotated,150)
assert inside(t,rotated), ('rotation',key)
counts['rotation'] += 1
evidence.append({'nu':str(nu),'x':str(x),'number':t,'integral_enclosure':str(integral)})
counts['rows'] += 1
if (index+1)%100 == 0:
print('checked',index+1,'of 600',flush=True)
bessel = {}
if stored:
bdoc = numberdb.table('T196')
bessel = {tuple(QQ(x) for x in identity.split(',')): written_ball(value)
for identity,value in _stored_entries(bdoc).items()}
for x in sorted({k[1] for k in keys}):
l0, l1 = values[(QQ(0),x)],values[(QQ(1),x)]
negative = struve_l(QQ(-1),x,600)
residual_ok(negative-2/R.pi()-l1,l1)
counts['negative_order'] += 1
if (QQ(0),x) in bessel and (QQ(1),x) in bessel:
i0,i1 = bessel[(QQ(0),x)],bessel[(QQ(1),x)]
sane(i0,95); sane(i1,95)
counts['stored_bessel_arguments'] += 1
else:
i0,i1 = real(C(x).bessel_I(0)),real(C(x).bessel_I(1))
sane(i0,150);sane(i1,150)
exact_integral,tail = integral_i0(x)
sane(exact_integral,150)
rhs = R.pi()*x/2*(i0*(2/R.pi()+l1)-i1*l0)
residual_ok(rhs-exact_integral,exact_integral)
maxima['i0_tail_bound'] = max(maxima['i0_tail_bound'],R(tail))
counts['bessel_integral'] += 1
if stored:
assert counts['stored_bessel_arguments'] == 30
summary = {'stored':stored, 'counts':counts,
'measurement':measure({str(k):v for k,v in text.items()}),
'maxima':{k:str(v) for k,v in maxima.items()},
'versions':{'numberdb':numberdb.__version__,'sage':sage.version.version,'mpmath':mp.__version__},
'working_bits':numberdb.bits(100,losing=64),
'quadrature_bits':480,'identity_bits':600}
print('EVIDENCE_JSON '+json.dumps(evidence),flush=True)
print('SUMMARY_JSON '+json.dumps(summary),flush=True)
return summary
if __name__ == '__main__':
if '--key-stdin' in sys.argv or os.environ.get('NUMBERDB_KEY_FROM_STDIN')=='1':
numberdb.configure(api_key=sys.stdin.read().strip())
main('--stored' in sys.argv)