back to table · edit · history · where entries came from · files · download
7198 bytes, as of the version from 2026-10-27 02:01 (current). Recorded here, not run.
"""Bessel moments, T521: the original quadrature extended to five even moments.
python -m pip install numberdb python-flint mpmath
python generate.py --compute --new-only
python generate.py --compute
python generate.py
The default enumeration contains all 40 original odd moments and five even
moments specified in even-selection.md. The last command verifies the live
table without writing. No command publishes, reviews, or approves anything.
This extends the attached generator with SHA256
6f50e2a959efd804b670232fe50fdb6ab25048e388278cc83652438a94aa9cfc:
the same direct tanh-sinh quadrature, Arb Bessel evaluations at exact binary
nodes, 150/200 working digits, radius gates, and exponential tail estimate.
Python-flint replaces Sage's Arb interface; no special-function method is
replaced. Old decimal strings are retained only after reproduction checks.
For k=0 only, t=u^2 changes the integrand to
2*u*I0(u^2)^a*K0(u^2)^b, whose limit at zero is zero. This avoids falsely
assigning a finite limiting value to K0(t)^b itself. Positive k uses the
original variable and panels. Both cases truncate at the same t=T.
For t>=T>=1, I0(t)<=exp(t), K0(t)<=sqrt(pi/(2*t))*exp(-t), b<=6,
b-a>=1 and 0<=k<=7 imply the common bound
(pi/2)^3 * 7! * exp(-T) * sum(T^j/j!, j=0..7)
for the omitted tail. These inequalities follow from DLMF 10.32.1 and
10.32.9, using cosh(v)>=1+v^2/2 for K0. T=3*(working+60).
The substitution does not change the integral or its tail at t=T.
Arb midpoints enter mpmath quadrature; discretization is NOT certified.
The honest label remains heuristic (agreement-checked), at 100 digits.
Independent gamma and elliptic-integral identities are in even-checks.py.
"""
import argparse
import json
from functools import lru_cache
import mpmath as mp
import numberdb
from flint import arb, ctx
WORKING_DIGITS = (150, 200)
STORED_DIGITS = 100
EVEN_KEYS = ((0, 3, 0), (0, 3, 2), (0, 4, 0), (0, 4, 2), (1, 3, 0))
def selection():
for total in range(3, 7):
for first in range((total - 1) // 2 + 1):
second = total - first
powers = {1, 3, 5, 7}
powers.update(power for left, right, power in EVEN_KEYS
if (left, right) == (first, second))
for power in sorted(powers):
yield {'a': first, 'b': second, 'k': power}
def binary_midpoint(value):
mantissa, exponent = value.man_exp()
return mp.ldexp(mp.mpf(int(mantissa)), int(exponent))
@lru_cache(maxsize=65536)
def _bessels(binary_node, working):
sign, mantissa, exponent, _ = binary_node
node = mp.mp.make_mpf(binary_node)
previous_precision = ctx.prec
ctx.prec = numberdb.bits(working + 40) + 3 * int(mp.ceil(node))
try:
exact_node = arb((-1 if sign else 1) * mantissa) * arb(2) ** exponent
values = []
for ball in (exact_node.bessel_i(0), exact_node.bessel_k(0)):
if not ball.is_finite():
raise ArithmeticError('Nonfinite Bessel enclosure')
midpoint = binary_midpoint(ball.mid())
radius = binary_midpoint(ball.rad())
if not mp.isfinite(midpoint) or midpoint <= 0:
raise ArithmeticError('Bessel midpoint is not finite and positive')
if not radius < midpoint * mp.mpf(10) ** (-working - 15):
raise ArithmeticError('Insufficient Bessel enclosure accuracy')
values.append(midpoint)
return tuple(values)
finally:
ctx.prec = previous_precision
def tail_bound(working):
with mp.workdps(working):
endpoint = 3 * (working + 60)
return ((mp.pi / 2) ** 3 * mp.factorial(7) * mp.exp(-endpoint)
* mp.fsum(mp.mpf(endpoint) ** power / mp.factorial(power)
for power in range(8)))
def integral(params, working):
first, second, power = (int(params[name]) for name in ('a', 'b', 'k'))
if not (0 <= first < second <= 6 and 0 <= power <= 7):
raise ValueError('This bound covers 0<=a<b<=6 and 0<=k<=7 only')
with mp.workdps(working):
def integrand(variable):
if not variable:
return mp.mpf(0)
node = variable ** 2 if power == 0 else variable
jacobian = 2 * variable if power == 0 else 1
first_value, second_value = _bessels(node._mpf_, working)
return jacobian * node ** power * first_value ** first * second_value ** second
endpoint = 3 * (working + 60)
panels = [0, 1, 4, 16, 64, endpoint]
if power == 0:
panels = [mp.sqrt(point) for point in panels]
value, estimated_error = mp.quad(integrand, panels, method='tanh-sinh', error=True)
if not mp.isfinite(value) or value <= 0:
raise ArithmeticError('Nonfinite or nonpositive moment')
if estimated_error >= abs(value) * mp.mpf(10) ** (-working + 10):
raise ArithmeticError('Quadrature convergence estimate failed')
if tail_bound(working) >= abs(value) * mp.mpf(10) ** (-working):
raise ArithmeticError('Tail bound too large')
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)
if relative >= mp.mpf('1e-130'):
raise ArithmeticError('150/200-digit agreement failed')
return {
'params': params, 'number': mp.nstr(high, digits, strip_zeros=False),
'at_150': mp.nstr(low, 150), 'at_200': mp.nstr(high, 200),
'relative_difference': mp.nstr(relative, 12),
'upper_tail_bounds': [mp.nstr(tail_bound(working), 12)
for working in WORKING_DIGITS],
'endpoint_substitution': 't=u^2' if int(params['k']) == 0 else 'none',
'rigour': 'heuristic (agreement-checked)'}
class BesselMoments(numberdb.Generator):
table = 'T521'
parameters = ('a', 'b', 'k')
type = 'R'
digits = STORED_DIGITS
rigour = 'heuristic (agreement-checked)'
files = ('generate.py', 'even-checks.py', 'even-selection.md')
def enumerate(self):
yield from selection()
def value(self, params, digits):
return evaluated(params, digits)['number']
if __name__ == '__main__':
parser = argparse.ArgumentParser(description=__doc__)
parser.add_argument('--compute', action='store_true')
parser.add_argument('--new-only', action='store_true')
options = parser.parse_args()
if options.compute:
for parameters in selection():
identity = tuple(parameters[name] for name in ('a', 'b', 'k'))
if not options.new_only or identity in EVEN_KEYS:
print(json.dumps(evaluated(parameters)), flush=True)
else:
if options.new_only:
parser.error('--new-only requires --compute')
report = BesselMoments().verify(sample=None)
print(report)
raise SystemExit(0 if report.ok else 1)