back to table · edit · history · where entries came from · files · download
8341 bytes, as of the version from 2026-10-26 14:52 (current). Recorded here, not run.
"""Independent checks for prolate differential eigenvalues -- numberdb.org/T518.
sage -python checks.py
sage -python checks.py --stored
The second command reads T518 and checks every stored decimal interval by
opposite signs of an independent endpoint-series residual. A private draft
requires NUMBERDB_API_KEY (or NUMBERDB_KEY_FROM_STDIN=1 and a piped key).
Output is JSON evidence, suitable for saving as evidence.json. No writes.
The endpoint solution is normalized here by f(1)=1, independently of the
unit-norm convention; eigenvalues do not depend on this normalization.
Writing t=1-x and f=sum a_k t^k in the differential equation gives
2(k+1)^2 a_(k+1) = (k(k+1)+q-chi)a_k - 2q a_(k-1) + q a_(k-2).
We impose f(0)=0 for odd modes or f'(0)=0 for even modes. This does not
construct a Legendre matrix or share its arithmetic. Fixed-precision series
and sign tests are diagnostic, not rigorous remainder bounds.
"""
import json
import os
import sys
import time
import platform
from pathlib import Path
import importlib.metadata
import numberdb.sage as numberdb
from sage.rings.rational_field import QQ
from sage.rings.real_mpfr import RealField
from scipy.special import pro_cv
import scipy
import sage.version
import generate as g
def endpoint(chi, q, parity, terms):
R = chi.parent()
a, previous, older = R(1), R(0), R(0)
da, dprevious, dolder = R(0), R(0), R(0)
total = R(1) if parity else R(0)
dtotal = R(0)
for k in range(terms-1):
factor = k*(k+1)+q-chi
divisor = 2*(k+1)**2
new = (factor*a - 2*q*previous + q*older)/divisor
dnew = (factor*da - a - 2*q*dprevious + q*dolder)/divisor
weight = 1 if parity else k+1
total += weight*new
dtotal += weight*dnew
older, previous, a = previous, a, new
dolder, dprevious, da = dprevious, da, dnew
return total, dtotal
def endpoint_root(S, n, dps, terms):
R = RealField(numberdb.bits(dps, losing=16))
q = (R.pi()*R(S)/2)**2
# Deliberately only a double-precision seed from an outside implementation.
chi = R(str(pro_cv(0, n, float(R.pi()*S/2))))
for iteration in range(12):
residual, derivative = endpoint(chi, q, n%2, terms)
assert derivative and (not derivative.is_infinity() and not derivative.is_NaN()) and (not residual.is_infinity() and not residual.is_NaN())
change = residual/derivative
chi -= change
if abs(change) < R(10)**(-(dps-30)):
return chi, iteration+1
raise ArithmeticError('Endpoint Newton iteration did not converge')
def agreement_digits(a, b):
R = RealField(900)
difference = abs(R(a)-R(b))
if not difference:
return 250
return int(-(difference/max(abs(R(a)),R(1))).log10())
def controls():
# Exact, nonzero controls before interpreting zero residuals.
assert endpoint(QQ(0), QQ(0), 1, 20)[0] == 1 # P0(0)
assert endpoint(QQ(2), QQ(0), 0, 20)[0] == -1 # d P1(1-t)/dt
assert endpoint(QQ(6), QQ(0), 1, 20)[0] == -QQ(1)/2 # P2(0)
for n in range(g.MODES):
residual, derivative = endpoint(QQ(n*(n+1)), QQ(0), n%2, 30)
assert residual == 0 and derivative != 0
computed = g.bisect_mode(g.field(160)(0), n, 20, 160)
assert abs(computed-n*(n+1)) < g.field(160)(10)**(-140)
# Algebraic check of the DLMF shift, in every degree used in the largest matrix.
for k in range(360):
diag = QQ((k+1)**2)/((2*k+1)*(2*k+3)) + QQ(k*k)/((2*k-1)*(2*k+1))
dlmf = -QQ(2*(k*(k+1)-1))/((2*k-1)*(2*k+3))
assert diag-1 == dlmf
upper = -QQ((k+1)*(k+2))/((2*k+3)*(2*k+5))
lower = -QQ((k+1)*(k+2))/((2*k+1)*(2*k+3))
square = QQ((k+1)**2*(k+2)**2)/((2*k+1)*(2*k+3)**2*(2*k+5))
assert upper*lower == square
return {'endpoint_exact_nonzero_controls':3,'zero_bandwidth_modes':12,
'dlmf_shift_rational_coefficient_checks':360}
def parse_decimal(text):
"""Exact rational midpoint and last decimal unit, without binary floats."""
mantissa, marker, exponent = text.lower().partition('e')
exponent = int(exponent) if marker else 0
integer, dot, fraction = mantissa.partition('.')
assert dot
power = exponent-len(fraction)
unit = QQ(10)**power
return QQ(int(integer+fraction, 10))*unit, unit
def stored_values():
tree = numberdb.table('T518')
values = {}
def walk(node, keys):
if len(keys)==3:
value = node['number'] if isinstance(node,dict) else node
values[tuple(map(int, keys))] = value
return
if isinstance(node,dict) and 'numbers' in node:
node=node['numbers']
for key, value in node.items():
walk(value, keys+[key])
nums = tree['Numbers']
if isinstance(nums,list):
for item in nums:
p=item['params']
values[(int(p['m']),int(p['S']),int(p['n']))]=item['number']
else:
walk(nums,[])
assert len(values)==48
return values
def run(use_stored=False):
start=time.monotonic()
report={'table':'T518','controls':controls(),'settings':g.SETTINGS,
'endpoint_settings':[(220,600),(260,800)],'rows':[],
'versions':{'python':platform.python_version(),'sage':sage.version.version,
'scipy':scipy.__version__}}
try:
report['versions']['numberdb']=importlib.metadata.version('numberdb')
except importlib.metadata.PackageNotFoundError:
report['versions']['numberdb']='repository client'
stored=stored_values() if use_stored else None
for S in g.BANDWIDTHS:
spectra=[g.spectrum(str(S),*setting) for setting in g.SETTINGS]
R=g.field(260)
q=(R.pi()*S/2)**2
previous=None
for n in range(g.MODES):
reference=R(spectra[-1][n])
settings_agreement=min(agreement_digits(v[n],reference) for v in spectra[:-1])
assert settings_agreement>=110
root1, it1=endpoint_root(S,n,220,600)
root2, it2=endpoint_root(S,n,260,800)
independent_agreement=min(agreement_digits(reference,root1),agreement_digits(reference,root2))
assert independent_agreement>=110
assert n*(n+1)<reference<n*(n+1)+q
if previous is not None:
assert reference>previous
previous=reference
scipy_error=abs(reference-R(str(pro_cv(0,n,float(R.pi()*S/2)))))/reference
assert scipy_error<R('1e-12')
text=stored[(0,S,n)] if stored is not None else g.decimal_string(reference,100)
midpoint,ulp=parse_decimal(text)
assert abs(R(midpoint)-root2)<R(ulp)
# Direct test of the interval asserted by the written decimal.
# At both series sizes, it brackets a root by opposite residual signs.
for dps,terms in ((220,600),(260,800)):
F=RealField(numberdb.bits(dps,losing=16))
qf=(F.pi()*S/2)**2
fl=endpoint(F(midpoint-ulp),qf,n%2,terms)[0]
fh=endpoint(F(midpoint+ulp),qf,n%2,terms)[0]
assert (not fl.is_infinity() and not fl.is_NaN()) and (not fh.is_infinity() and not fh.is_NaN()) and fl*fh<0
report['rows'].append({'m':0,'S':S,'n':n,'number':text,
'precision_truncation_agree_digits':settings_agreement,
'endpoint_agree_digits':independent_agreement,
'endpoint_self_agree_digits':agreement_digits(root1,root2),
'endpoint_iterations':[it1,it2],
'scipy_relative_error':str(scipy_error.n(30)),
'decimal_interval_sign_changes':2,
'rounding_error_ulps':str((abs(R(midpoint)-root2)/R(ulp)).n(30))})
print('checked S=%d'%S,file=sys.stderr,flush=True)
report['read_from_database']=use_stored
report['minimum_precision_truncation_agreement']=min(r['precision_truncation_agree_digits'] for r in report['rows'])
report['minimum_endpoint_agreement']=min(r['endpoint_agree_digits'] for r in report['rows'])
report['seconds']=time.monotonic()-start
return report
if __name__=='__main__':
if os.environ.get('NUMBERDB_KEY_FROM_STDIN')=='1':
os.environ['NUMBERDB_API_KEY']=sys.stdin.read().strip()
print(json.dumps(run('--stored' in sys.argv),indent=2))