back to table · edit · history · where entries came from · files · download
4625 bytes, as of the version from 2026-10-26 16:56. Recorded here, not run.
"""Bessel moments -- numberdb.org/T521.
Run with SageMath (or passagemath with numberdb and mpmath):
sage -pip install numberdb mpmath
sage -python generate.py # verify all stored entries
sage -python generate.py --compute # print reproducible computations
The default selection is explained in argument-selection.md. Direct tanh-sinh
quadrature at 150 and 200 decimal working digits produces 100 significant
stored digits. Arb evaluates the Bessel functions, but its midpoints enter
mpmath quadrature: the claim is heuristic (agreement-checked), not certified
quadrature. The same nodes' Bessel evaluations are cached, never the moments.
No recurrence or closed-form identity generates the entries.
The upper cutoff is T=3*(working+60). For these 40 integrals, t>=T>=1,
I0(t)<=exp(t) and K0(t)<=sqrt(pi/(2*t))*exp(-t) imply a common tail bound
(pi/2)^3 * 7! * exp(-T) * sum(T^j/j!, j=0..7).
The lower endpoint is integrated by tanh-sinh using the limiting value zero;
t^k*log(t)^b tends to zero for the selected positive k. Quadrature convergence
and agreement do not certify its discretization error.
"""
import numberdb.sage as numberdb
from sage.rings.complex_arb import ComplexBallField
from sage.rings.rational_field import QQ
from functools import lru_cache
import json
import mpmath as mp
import sys
WORKING_DIGITS = (150, 200)
STORED_DIGITS = 100
def selection():
for total in range(3, 7):
for a in range((total - 1) // 2 + 1):
for k in (1, 3, 5, 7):
yield {'a': a, 'b': total - a, 'k': k}
@lru_cache(maxsize=65536)
def _bessels(binary_x, working):
# Preserve the exact binary quadrature node when passing it to Arb.
sign, mantissa, exponent, _ = binary_x
q = QQ((-1 if sign else 1) * mantissa) * QQ(2)**exponent
# Allow for exp(t) versus exp(-t) cancellation in the Bessel evaluator.
# 3*t bits exceeds 2*t/log(2); each returned radius is also checked.
field = ComplexBallField(numberdb.bits(working + 40) + 3 * int(q.ceil()))
x = field(q)
iv, kv = x.bessel_I(0), x.bessel_K(0)
assert iv.imag().contains_zero() and kv.imag().contains_zero()
values = []
for v in (iv.real(), kv.real()):
assert v.is_finite()
m = mp.mpf(str(v.mid()))
assert mp.isfinite(m) and m > 0
assert mp.mpf(str(v.rad())) < m * mp.mpf(10)**(-working - 15), (working, str(q)[:80], mp.nstr(m, 15), str(v.rad()))
values.append(m)
return tuple(values)
def integral(params, working):
a, b, k = (int(params[name]) for name in ('a', 'b', 'k'))
if not (0 <= a < b and k >= 1):
raise ValueError('This quadrature implements a>=0, b>a, k>=1.')
with mp.workdps(working):
def f(t):
if not t:
return mp.mpf(0)
iv, kv = _bessels(t._mpf_, working)
return t**k * iv**a * kv**b
endpoint = 3 * (working + 60)
value, estimated_error = mp.quad(f, [0, 1, 4, 16, 64, endpoint],
method='tanh-sinh', error=True)
assert mp.isfinite(value) and value > 0
assert estimated_error < abs(value) * mp.mpf(10)**(-working + 10)
return +value
def evaluated(params, digits=STORED_DIGITS):
if digits > STORED_DIGITS:
raise ValueError('This run justifies at most 100 stored digits.')
low = integral(params, WORKING_DIGITS[0])
high = integral(params, WORKING_DIGITS[1])
with mp.workdps(WORKING_DIGITS[1]):
relative = abs(low - high) / abs(high)
assert relative < mp.mpf(10)**(-130), (params, relative)
text = mp.nstr(high, digits, strip_zeros=False)
return text, mp.nstr(low, 150), mp.nstr(high, 200), mp.nstr(relative, 12)
class BesselMoments(numberdb.Generator):
table = 'T521'
parameters = ('a', 'b', 'k')
type = 'R'
digits = STORED_DIGITS
rigour = 'heuristic (agreement-checked)'
files = ('generate.py', 'checks.py', 'argument-selection.md')
def enumerate(self):
yield from selection()
def value(self, params, digits):
return evaluated(params, digits)[0]
if __name__ == '__main__':
if '--compute' in sys.argv:
for params in selection():
value, low, high, relative = evaluated(params)
print(json.dumps({'params': params, 'number': value,
'at_150': low, 'at_200': high,
'relative_difference': relative}), flush=True)
else:
report = BesselMoments().verify(sample=None)
print(report)
sys.exit(0 if report.ok else 1)