back to table · edit · history · where entries came from · files · download
8239 bytes, as of the version from 2026-10-27 01:18 (current). Recorded here, not run.
"""Independent build checks for Bessel Rayleigh sums -- numberdb.org/T533.
Run: sage -python checks.py
Read the stored draft and check it: sage -python checks.py --stored
The private draft needs NUMBERDB_API_KEY or --key-stdin. No calls write.
Generation uses normalized Bessel Taylor coefficients. This file separately
uses Kishore's convolution recurrence (Gupta and Muldoon, equations (3)-(4)),
Bernoulli numbers (sine-spectrum specialization), and mpmath Bessel roots.
Root sums are supplementary numerical checks, not exactness certificates.
Their omitted tails use j_(nu,k) > (k-1/2)*pi for nu>=0, by order monotonicity
and J_(-1/2)(z)=sqrt(2/(pi*z))*cos(z), then the decreasing integral test.
The roots are computed at 80 decimal digits, not certified as indexed roots.
Instruction evidence: agents/runs/20261026T235354Z-build-skills/manifest.json;
both archived skills have SHA256
b4d51ded0d52cfe3dcc0859de0845c3879b0c2363273ba7e809a7e06ded40b18.
"""
import numberdb.sage as numberdb
from sage.rings.rational_field import QQ
from sage.arith.misc import bernoulli
from fractions import Fraction
from math import factorial
from functools import lru_cache
from pathlib import Path
import json
import sys
import re
import hashlib
import platform
import mpmath as mp
from generate import BesselRayleighSums
from check import exactness, measure, agrees_with, names_its_rings
@lru_cache(maxsize=None)
def convolution(order, maximum):
"""Separate implementation and arithmetic library from generation."""
a = Fraction(order)
s = [Fraction(0), Fraction(1, 4) / (a + 1)]
for n in range(2, maximum + 1):
s.append(sum((s[k] * s[n-k] for k in range(1, n)), Fraction(0)) / (a+n))
return tuple(s[1:])
def rational(text):
f = Fraction(str(text))
return QQ(f.numerator) / QQ(f.denominator)
def other(key):
order, n = key
return rational(convolution(order, n)[-1])
def flatten_stored(document):
block = document['Numbers']
if isinstance(block, list):
return {(str(r['params']['nu']), int(r['params']['n'])): rational(r['number'])
for r in block}
return {(str(nu), int(n)): rational(v['number'] if isinstance(v, dict) else v)
for nu, group in block.items() for n, v in group.items()}
def numerical_roots(values):
mp.mp.dps = 80
count = 32
tolerance = mp.mpf('1e-65')
# Nonzero known controls precede the measurements.
for k in range(1, count + 1):
z = mp.besseljzero(mp.mpf('0.5'), k)
assert abs(z - k*mp.pi) < tolerance
control = sum((1/(k*mp.pi)**2 for k in range(1, count+1)), mp.mpf(0))
control_tail = 1/(mp.pi**2 * (count-mp.mpf('0.5')))
assert control > mp.mpf('0.16')
assert control < mp.mpf(1)/6 < control + control_tail
def in_tail_interval(target, partial, tail):
return tolerance < target-partial < tail-tolerance
rows = []
orders = list(dict.fromkeys(order for order, _ in values))
for order in orders:
nu = mp.mpf(Fraction(order).numerator)/Fraction(order).denominator
roots = [mp.besseljzero(nu, k) for k in range(1, count+1)]
assert all(roots[k] < roots[k+1] for k in range(count-1))
for n in range(1, 13):
q = values[(order, n)]
target = mp.mpf(str(q.numerator())) / mp.mpf(str(q.denominator()))
partial = sum((z**(-2*n) for z in roots), mp.mpf(0))
tail = (mp.pi**(-2*n) * (count-mp.mpf('0.5'))**(1-2*n)
/ (2*n-1))
gap = target-partial
assert mp.isfinite(partial) and mp.isfinite(tail)
assert in_tail_interval(target, partial, tail), (order, n, gap, tail)
rows.append({'nu': order, 'n': n, 'partial_sum': mp.nstr(partial, 70),
'tail_bound': mp.nstr(tail, 15),
'gap_over_bound': mp.nstr(gap/tail, 15),
'relative_tail_bound': mp.nstr(tail/target, 15)})
# The test must reject a bad value; adding twice the tail leaves the bracket.
assert not in_tail_interval(target + 2*tail, partial, tail)
return {'root_count_per_order':count,'working_decimal_digits':80,
'control_half_order_roots':32,'known_sum_control':'1/6',
'checked_entries':len(rows),'rows':rows,
'limitation':'Numerically indexed roots; this check does not certify 100 digits. Exactness comes from rational coefficient extraction and the canonical product.'}
def run(values, with_roots=True):
expected = [(a, n) for a in ('0','1/2','1','3/2','2','5/2','3','7/2')
for n in range(1, 13)]
assert list(values) == expected
assert not exactness(values)
assert not names_its_rings(Path(__file__).with_name('generate.py'))
# Known nonzero controls for the separate rational method.
assert other(('0',1)) == QQ(1)/QQ(4)
assert other(('1/2',2)) == QQ(1)/QQ(90)
assert not agrees_with(values, other)
perturbed = dict(values)
perturbed[('7/2',12)] += QQ(1)/QQ(10**80)
assert len(agrees_with(perturbed, other)) == 1
low_count = 0
for a in ('0','1/2','1','3/2','2','5/2','3','7/2'):
v = rational(a)
formulas = [QQ(1)/(4*(v+1)),
QQ(1)/(16*(v+1)**2*(v+2)),
QQ(1)/(32*(v+1)**3*(v+2)*(v+3))]
for n, answer in enumerate(formulas,1):
assert values[(a,n)] == answer
low_count += 1
for n in range(1,13):
assert values[(a,n)] > 0
for n in range(1,12):
assert values[(a,n)] > values[(a,n+1)]
for n in range(1,11):
assert values[(a,n)]*values[(a,n+2)] > values[(a,n+1)]**2
assert bernoulli(2) == QQ(1)/QQ(6)
assert bernoulli(4) == -QQ(1)/QQ(30)
for n in range(1,13):
expected_half = (QQ((-1)**(n+1) * 2**(2*n-1)) * bernoulli(2*n)
/ QQ(factorial(2*n)))
assert values[('1/2',n)] == expected_half
# Independent transcendental evaluation of the same specialization.
with mp.workdps(120):
zeta_ratio = mp.zeta(2*n)/mp.pi**(2*n)
rat = values[('1/2',n)]
exact_as_mp = mp.mpf(str(rat.numerator()))/mp.mpf(str(rat.denominator()))
assert abs(zeta_ratio/exact_as_mp-1) < mp.mpf('1e-115')
result = {'entries':len(values),'convolution_matches':len(values),
'low_order_formula_matches':low_count,'bernoulli_specializations':12,
'zeta_specializations_120_dps':12,'positivity':96,
'strict_decrease':88,'strict_log_convexity':80,
'mutation_control':'detected one changed entry at (7/2,12)',
'measure':measure(values)}
if with_roots:
result['root_checks'] = numerical_roots(values)
return result
if __name__ == '__main__':
client = None
if '--key-stdin' in sys.argv:
key = sys.stdin.read().strip()
if re.match(r'^[A-Z_][A-Z0-9_]*\s*=',key):
key = key.split('=',1)[1].strip().strip('"\'')
client = numberdb.Client(api_key=key)
generator = BesselRayleighSums()
if '--stored' in sys.argv:
values = flatten_stored(numberdb.table(generator.table,client=client))
report = generator.verify(sample=None,client=client)
assert report.ok, report
print(report, file=sys.stderr)
else:
values = {(p['nu'],int(p['n'])):generator.value(p,100)
for p in generator.enumerate()}
evidence = run(values)
evidence['mode'] = 'stored' if '--stored' in sys.argv else 'generated'
evidence['skill_manifest'] = 'agents/runs/20261026T235354Z-build-skills/manifest.json'
evidence['skill_sha256'] = 'b4d51ded0d52cfe3dcc0859de0845c3879b0c2363273ba7e809a7e06ded40b18'
evidence['versions'] = {'python':platform.python_version(),'mpmath':mp.__version__,
'generator_environment':generator.environment()}
evidence['sha256'] = {name:hashlib.sha256(Path(__file__).with_name(name).read_bytes()).hexdigest()
for name in ['generate.py','checks.py','check.py']}
evidence['values'] = [{'nu':a,'n':n,'number':str(v)} for (a,n),v in values.items()]
print(json.dumps(evidence,indent=2))