back to table · edit · history · where entries came from · files · download
8542 bytes, as of the version from 2026-10-27 03:18 (current). Recorded here, not run.
"""Independent checks for Rayleigh sums of Bessel derivatives -- numberdb.org/T535.
sage -python checks.py checks the generated selection
sage -python checks.py --stored checks database values, needs draft API access
Only --stored reads the website; neither command writes. The exact reference
computation uses Python Fraction and the Riccati recurrence obtained from
Gupta-Muldoon (math/9910128), equations (19)-(21), a=c=0, b=1. The upper limits
in their specialized equation (23) omit endpoint terms; we do not use it.
The generator instead uses a linear Newton recurrence on Taylor coefficients.
Numerical root sums are supplemental checks, not certificates of the exact
values or 100-digit enclosures. Interlacing and monotonicity (DLMF 10.21(i),(iv))
and j_(1/2,k)=k*pi imply j'_(nu,k)>(k-1)*pi for nu=1,2,3, k>=2.
For nu=0 the positive roots equal j_(1,k)>k*pi. Thus for K>=2, the omitted
tail is at most pi^(-2*n)*(K-1)^(1-2*n)/(2*n-1). Approximate root indexing
and rounding in mpmath are not rigorous; their error is tested by controls,
residuals and a precision rerun rather than asserted to be an enclosure.
"""
import json
import os
import sys
from fractions import Fraction as F
from pathlib import Path
import platform
import numberdb.sage as numberdb
from sage.rings.rational_field import QQ
import sage.version
import mpmath as mp
from generate import BesselDerivativeRayleigh
from check import exactness, measure, agrees_with, names_its_rings
def bessel_moments(nu, limit):
nu = F(nu)
s = [F(0), F(1,4)/(nu+1)]
for n in range(2, limit+1):
s.append(sum(s[k]*s[n-k] for k in range(1,n))/(nu+n))
return s
def riccati(nu, limit):
if nu == 0:
return bessel_moments(1, limit)
v = F(nu)
s = [F(0), (v+2)/(4*v*(v+1))]
if limit > 1:
s.append((v*v+8*v+8)/(16*v*v*(v+1)**2*(v+2)))
for n in range(3,limit+1):
quadratic = sum(s[k]*s[n-k] for k in range(1,n))
shifted = sum(s[k]*s[n-1-k] for k in range(1,n-1))
s.append(((n+v-2)*s[n-1]+v*v*quadratic-shifted)/(v*v*(n+v)))
return s
def expected_values():
return {(nu,n):riccati(nu,12)[n] for nu in range(4) for n in range(1,13)}
def read_values():
g = BesselDerivativeRayleigh()
if '--stored' not in sys.argv:
raw = {(int(p['nu']),int(p['n'])):g.value(p,100) for p in g.enumerate()}
assert not exactness(raw)
# Compare the two arithmetic systems through exact integer pairs.
return {key:F(int(v.numerator()),int(v.denominator())) for key,v in raw.items()}
if os.environ.get('NUMBERDB_KEY_FROM_STDIN') == '1':
key = sys.stdin.read().strip()
if '=' in key and key.split('=',1)[0].replace('_','').isalnum():
key = key.split('=',1)[1].strip().strip('\"\'')
numberdb.configure(api_key=key)
report = g.verify(sample=None)
assert report.ok, str(report)
document = numberdb.table(g.table)
numbers = document['Numbers']
result = {}
if isinstance(numbers,list):
for row in numbers:
p = row['params']
result[(int(p['nu']),int(p['n']))] = F(row['number'])
else:
for nu, group in numbers.items():
for n, value in group.items():
result[(int(nu),int(n))] = F(value['number'] if isinstance(value,dict) else value)
print('verify(sample=None): all 48 stored entries agree', flush=True)
return result
def test_exact(values):
assert set(values) == {(v,n) for v in range(4) for n in range(1,13)}
assert all(isinstance(v,F) for v in values.values())
assert not exactness({k:QQ(v.numerator)/QQ(v.denominator) for k,v in values.items()})
assert not names_its_rings(Path(__file__).with_name('generate.py'))
expected = expected_values()
# Nonzero controls before trusting the independent recurrence.
assert riccati(0,3)[1:] == [F(1,8),F(1,192),F(1,3072)]
assert riccati(1,3)[1:] == [F(3,8),F(17,192),F(79,3072)]
assert bessel_moments(F(1,2),2)[1:] == [F(1,6),F(1,90)]
assert not agrees_with(values,lambda key: expected[key])
bad = dict(values)
bad[(1,12)] += F(1,10**30)
assert len(agrees_with(bad,lambda key: expected[key])) == 1
low = 0
for nu in (1,2,3):
v = F(nu)
forms = [(v+2)/(4*v*(v+1)),
(v*v+8*v+8)/(16*v*v*(v+1)**2*(v+2)),
(v**3+16*v*v+38*v+24)/(32*v**3*(v+1)**3*(v+2)*(v+3))]
for n, value in enumerate(forms,1):
assert values[(nu,n)] == value
low += 1
sigma = bessel_moments(nu,12)
for n in range(2,13):
assert values[(nu,n)] == -sigma[n]+F(2,nu)*sum(
sigma[k]*values[(nu,n-k)] for k in range(1,n))
for nu in range(4):
for n in range(1,13): assert values[(nu,n)] > 0
for n in range(1,12): assert values[(nu,n)] > values[(nu,n+1)]
for n in range(2,12):
assert values[(nu,n)]**2 < values[(nu,n-1)]*values[(nu,n+1)]
return {'independent_exact_matches':48, 'closed_forms_positive_order':low,
'mixed_Bessel_recurrence_matches':33,'known_controls':8,
'deliberate_perturbation_rejected':True,'positive':48,
'decreasing':44,'strict_log_convexity':40}
def numerical_check(values, precision):
mp.mp.dps = precision
K = 32
slack = mp.mpf(10)**(-precision+15)
control_roots = [mp.besseljzero(mp.mpf('0.5'),k) for k in range(1,K+1)]
assert max(abs(x-k*mp.pi) for k,x in enumerate(control_roots,1)) < slack
control_sum = mp.fsum(x**-2 for x in control_roots)
assert 0 < mp.mpf(1)/6-control_sum < 1/(mp.pi**2*(K-1))
assert abs(mp.besseljzero(0,1,derivative=1)) == 0
root_rows, comparison_rows = [], []
max_residual = mp.mpf(0)
for nu in range(4):
roots = [mp.besseljzero(nu,k+(1 if nu==0 else 0),derivative=1)
for k in range(1,K+1)]
assert all(x>0 for x in roots)
assert all(a<b for a,b in zip(roots,roots[1:]))
for k,x in enumerate(roots,1):
residual = abs(mp.besselj(nu,x,derivative=1))
max_residual = max(max_residual,residual)
assert residual < slack
assert x > (k-1)*mp.pi
root_rows.append({'nu':nu,'k':k,'root':mp.nstr(x,precision-10)})
for n in range(1,13):
partial = mp.fsum(x**(-2*n) for x in roots)
tail = (K-1)**(1-2*mp.mpf(n))/((2*n-1)*mp.pi**(2*n))
v = F(str(values[(nu,n)]))
exact = mp.mpf(v.numerator)/v.denominator
remainder = exact-partial
assert -slack < remainder < tail+slack
if nu==1 and n==1:
assert not (-slack < exact+mp.mpf('0.01')-partial < tail+slack)
comparison_rows.append({'nu':nu,'n':n,'partial_sum':mp.nstr(partial,precision-10),
'remainder':mp.nstr(remainder,12),'tail_bound':mp.nstr(tail,12)})
return {'dps':precision,'roots':root_rows,'comparisons':comparison_rows,
'max_residual':mp.nstr(max_residual,12),'root_indexing_certified':False,
'known_half_order_control':True,'zero_mode_control':True,
'deliberate_root_sum_perturbation_rejected':True}
def main():
values = read_values()
exact_report = test_exact(values)
print(json.dumps(exact_report), flush=True)
numerical = [numerical_check(values,p) for p in (70,100)]
for a,b in zip(numerical[0]['roots'],numerical[1]['roots']):
assert abs(mp.mpf(a['root'])-mp.mpf(b['root'])) < mp.mpf('1e-56')
if '--stored' in sys.argv:
other = numberdb.table('T533')['Numbers']['1']
assert all(values[(0,n)] == F(other[str(n)]) for n in range(1,13))
exact_report['order_zero_matches_T533_order_one'] = 12
result = {'table':'T535','mode':'stored' if '--stored' in sys.argv else 'local',
'skill_manifest':'agents/runs/20261027T024416Z-build-skills/manifest.json',
'skill_sha256':'e37efb7e2ededf6c1d7134c673f152f648f26e0267c556aebdf3ec65765d8b4c',
'versions':{'python':platform.python_version(),'sage':sage.version.version,
'numberdb':numberdb.__version__,'mpmath':mp.__version__},
'exact_checks':exact_report,'measurement':measure(values),
'values':[{'nu':nu,'n':n,'number':str(value)} for (nu,n),value in values.items()],
'numerical':numerical}
print('EVIDENCE_JSON='+json.dumps(result,default=str), flush=True)
if __name__ == '__main__':
main()