back to table · edit · history · where entries came from · files · download
6503 bytes, as of the version from 2026-10-27 02:01 (current). Recorded here, not run.
"""Independent numerical checks for Bessel moments -- numberdb.org/T521.
Run: sage -python checks.py computation.jsonl
Or: sage -python checks.py stored-document.json
No generator functions are imported. Published recurrences, zeta/L-values,
and exact multinomial counts check the digits. A second quadrature uses
mpmath's own Bessel implementation, Gauss-Legendre nodes and t=exp(x), rather
than the generator's Arb Bessel values and direct tanh-sinh quadrature.
The nonzero control IKM(0,2;1)=1/2 is checked first. All 40 rows are compared.
The independent quadrature uses 60 working digits and a 50-digit tolerance;
it does not claim an enclosure or independently certify all 100 stored digits.
"""
import numberdb.sage as numberdb
from functools import lru_cache
from math import factorial, comb
from pathlib import Path
import json
import mpmath as mp
import sys
import time
PAIRS = ((0,3),(1,2),(0,4),(1,3),(0,5),(1,4),(2,3),(0,6),(1,5),(2,4))
KEYS = [(a,b,k) for a,b in PAIRS for k in (1,3,5,7)]
def load_values(path):
raw = Path(path).read_text()
if raw.lstrip().startswith('{'):
try:
doc = json.loads(raw)
except json.JSONDecodeError:
doc = None
if doc is not None and 'Numbers' in doc:
result = {}
block = doc['Numbers']
if isinstance(block, list):
for entry in block:
key = tuple(int(entry['params'][n]) for n in ('a','b','k'))
result[key] = entry['number']
elif doc.get('Display properties', {}).get('group parameters') == [['a','b'],['k']]:
for pair, ks in block.items():
a,b = (int(x.strip()) for x in pair.split(','))
for k, entry in ks.items():
result[(a,b,int(k))] = entry['number'] if isinstance(entry,dict) else entry
else:
for a, bs in block.items():
for b, ks in bs.items():
for k, entry in ks.items():
result[(int(a),int(b),int(k))] = entry['number'] if isinstance(entry,dict) else entry
return result, False
rows = [json.loads(line) for line in raw.splitlines() if line.startswith('{')]
return {tuple(row['params'][n] for n in ('a','b','k')):
row.get('at_200',row['number']) for row in rows if 'params' in row}, True
def multinomial_square_sum(r, parts):
def compositions(total, length):
if length == 1:
yield (total,)
else:
for j in range(total+1):
for rest in compositions(total-j,length-1):
yield (j,)+rest
total=0
for composition in compositions(r,parts):
denominator=1
for j in composition:
denominator*=factorial(j)
coefficient=factorial(r)//denominator
total+=coefficient*coefficient
return total
def identities(values, precise):
assert set(values)==set(KEYS), 'Unexpected parameter selection'
with mp.workdps(220):
v={key:mp.mpf(val) for key,val in values.items()}
tolerance=mp.mpf('1e-180' if precise else '2e-98')
records=[]
def compare(label,terms):
residual=abs(mp.fsum(terms))/max(abs(t) for t in terms)
assert mp.isfinite(residual) and residual<tolerance,(label,mp.nstr(residual))
records.append({'check':label,'relative_residual':mp.nstr(residual,12)})
for a,b in PAIRS:
total=a+b
for k in ((1,3) if total in (3,4) else (1,)):
if total==3:
c=[(k+1)**4,-2*(5*k*k+20*k+21),9]
elif total==4:
c=[(k+1)**5,-4*(k+2)*(5*k*k+20*k+23),64*(k+3)]
elif total==5:
c=[(k+1)**6,-(35*(k+2)**4+42*(k+2)**2+3),259*(k+3)**2+104,-225]
else:
c=[(k+1)**7,-(k+2)*(56*(k+2)**4+112*(k+2)**2+24),
(k+3)*(784*(k+3)**2+944),-2304*(k+4)]
compare(f'recurrence {a},{b};{k}',[cj*v[a,b,k+2*j] for j,cj in enumerate(c)])
l2=(mp.zeta(2,mp.mpf(1)/3)-mp.zeta(2,mp.mpf(2)/3))/9
compare('IKM(0,3;1)=3L(2,chi_-3)/4',[v[0,3,1],-3*l2/4])
compare('IKM(0,4;1)=7zeta(3)/8',[v[0,4,1],-7*mp.zeta(3)/8])
compare('IKM(0,4;3)=7zeta(3)/32-3/16',[v[0,4,3],-7*mp.zeta(3)/32,mp.mpf(3)/16])
assert [multinomial_square_sum(r,3) for r in range(4)]==[1,3,15,93]
assert [multinomial_square_sum(r,4) for r in range(4)]==[1,4,28,256]
for r in range(4):
s3=mp.pi/(3*mp.sqrt(3))*(mp.mpf(2)**r*factorial(r)/mp.mpf(3)**r)**2*multinomial_square_sum(r,3)
s4=mp.pi**2/16*(mp.mpf(factorial(r))/4**r)**2*multinomial_square_sum(r,4)
compare(f'3-factor multinomial r={r}',[v[1,2,2*r+1],-s3])
compare(f'4-factor multinomial r={r}',[v[1,3,2*r+1],-s4])
assert all(x>0 and mp.isfinite(x) for x in v.values())
return records
@lru_cache(maxsize=32768)
def alternative_bessels(binary_x):
x=mp.mp.make_mpf(binary_x)
t=mp.exp(x)
return t,mp.besseli(0,t),mp.besselk(0,t)
def alternative(params):
a,b,k=params
with mp.workdps(60):
def f(x):
t,iv,kv=alternative_bessels(x._mpf_)
return t**(k+1)*iv**a*kv**b
return mp.quad(f,[-100,-16,-4,0,2,6],method='gauss-legendre')
def quadrature_checks(values):
start=time.time()
control=alternative((0,2,1))
with mp.workdps(70):
error=abs(control-mp.mpf(1)/2)
assert error<mp.mpf('1e-55'),error
print(json.dumps({'alternative_control':'IKM(0,2;1)=1/2','absolute_error':mp.nstr(error,12)}),flush=True)
for key in KEYS:
val=alternative(key)
with mp.workdps(70):
wanted=mp.mpf(values[key]);error=abs(val-wanted)/abs(wanted)
assert mp.isfinite(val) and error<mp.mpf('1e-50'),(key,error)
print(json.dumps({'alternative_params':key,'number':mp.nstr(val,60),
'relative_difference':mp.nstr(error,12),
'elapsed_seconds':round(time.time()-start,2)}),flush=True)
if __name__=='__main__':
values,precise=load_values(sys.argv[1])
records=identities(values,precise)
print(json.dumps({'identities':records,'count':len(records)}),flush=True)
quadrature_checks(values)
print(json.dumps({'independent_checks_complete':True,'entries':len(values)}),flush=True)