"""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 (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()