back to table · edit · history · where entries came from · files · download
3367 bytes, as of the version from 2026-10-26 15:32 (current). Recorded here, not run.
"""Checks for the presentation repair of numberdb.org/T519.
Run with SageMath, numberdb, NumPy and SciPy:
sage -python repair_checks.py T519-document.json
No API writes. The optional argument is an API document saved as JSON.
Checks kernel commutation and the trace exactly, plus independent double-
precision Nystrom checks. These numerical checks do not certify 100 digits.
"""
import json
import platform
import sys
from pathlib import Path
import numberdb.sage as numberdb
import sympy as sp
from sage.env import SAGE_VERSION
import mpmath as mp
import numpy as np
import scipy
from scipy.special import roots_legendre
from scipy.linalg import eigvalsh
x,y=sp.symbols('x y', real=True)
c=sp.symbols('c', positive=True)
kernel=sp.sin(c*(x-y))/(sp.pi*(x-y))
def differential(expr,z):
return -((1-z*z)*expr.diff(z)).diff(z)+c*c*z*z*expr
commutator=sp.simplify(differential(kernel,x)-differential(kernel,y))
assert commutator==0
assert (1-y*y).subs(y,1)==0 and (1-y*y).subs(y,-1)==0
assert sp.simplify(sp.limit(kernel,x,y)-c/sp.pi)==0
assert sp.simplify(sp.integrate(c/sp.pi,(x,-1,1))-2*c/sp.pi)==0
path=Path(sys.argv[1]) if len(sys.argv)>1 else Path(__file__).with_name('T519-document.json')
doc=json.loads(path.read_text()); mp.mp.dps=150
nodes,weights=roots_legendre(160)
constant=np.sqrt(weights[:,None]*weights[None,:])
control=eigvalsh(constant)
assert abs(control[-1]-2)<1e-12 and np.max(np.abs(control[:-1]))<1e-12
moment_error=max(abs(np.dot(weights,nodes**k)-(2/(k+1) if k%2==0 else 0)) for k in range(10))
assert moment_error<1e-12
rows=[]; summaries={}; largest_error=0
for S,modes in doc['Numbers'].items():
stored=[mp.mpf(modes[str(n)]) for n in range(12)]
assert all(0<v<1 for v in stored)
assert all(a>b for a,b in zip(stored,stored[1:]))
assert mp.fsum(stored)<int(S)
half=int(S)/2
operator=constant*half*np.sinc(half*(nodes[:,None]-nodes[None,:]))
computed=eigvalsh(operator)[::-1]
compared=0
for n,v in enumerate(stored):
if v>mp.mpf('1e-7'):
error=abs(computed[n]-float(v)); largest_error=max(largest_error,error)
assert error<1e-12, (S,n,error)
rows.append({'S':int(S),'n':n,'absolute_error':error}); compared+=1
summaries[S]={'ordered_and_in_unit_interval':True,'partial_trace_deficit':mp.nstr(int(S)-mp.fsum(stored),30),'modes_compared':compared}
exec(doc['Programs']['Python']['code'],ns:={})
program_value=float(ns['largest_concentration'](2))
assert abs(program_value-float(doc['Numbers']['2']['0']))<1e-12
print(json.dumps({'table':'T519','versions':{'python':platform.python_version(),'sage':SAGE_VERSION,'sympy':sp.__version__,'numberdb':numberdb.__version__,'numpy':np.__version__,'scipy':scipy.__version__,'mpmath':mp.__version__},'exact_checks':{'kernel_commutator':str(commutator),'endpoint_boundary_factors_zero':True,'diagonal_limit':'c/pi','trace':'2*c/pi'},'constant_kernel_max_error':float(max(abs(control[-1]-2),np.max(np.abs(control[:-1])))),'moment_max_error':float(moment_error),'all_stored_values_consistency_checked':48,'independent_values_compared':len(rows),'maximum_absolute_error':float(largest_error),'rows':rows,'bandwidth_checks':summaries,'displayed_program_value':program_value,'precision_limitation':'Independent comparisons use double precision; final stored digits are not re-certified by this repair.'},indent=2))