back to table · edit · history · where entries came from · files · download
12072 bytes, as of the version from 2026-10-27 02:15 (current). Recorded here, not run.
"""Independent numerical checks for PSWF zeros -- numberdb.org/T534.
sage -python checks.py > evidence.json
sage -python checks.py --stored > stored-evidence.json
Keep generate.py, spectral.py, argument-selection.md and skill-manifest.json
beside this file. --stored requires authenticated access for a private draft.
All operations are read-only. Each output records its sources and settings.
This checker solves the ODE from the endpoint x=1 with f(1)=1. Its Taylor
recurrence in t=1-x shares no Legendre matrix or eigenvector calculation with
the generator. It obtains chi by imposing even/odd parity at x=0, starting
from SciPy pro_cv; it then finds every positive root from this independent
function. Taylor lengths and working precision are both varied. These checks
have no rigorous infinite-series remainder or root-index certificate.
The endpoint recurrence and its chi derivative are adapted from T518's checks,
https://numberdb.org/files/T518/checks.py?raw=1 (downloaded 2026-10-27).
Root evaluation, interval tests, parity and zero-count checks are added here.
"""
import os
import sys
import json
import time
import hashlib
import platform
from pathlib import Path
from functools import lru_cache
import importlib.metadata
import numberdb.sage as numberdb
from sage.rings.rational_field import QQ
from sage.rings.real_mpfr import RealField
import sage.version
import scipy
from scipy.special import pro_cv, pro_ang1
import generate as g
ENDPOINT_SETTINGS=((220,600),(260,800))
def endpoint_coefficients(chi,q,terms,with_derivative=False):
R=chi.parent()
coeff=[R(1)]
derivative=[R(0)]
for k in range(terms-1):
a=coeff[k]
previous=coeff[k-1] if k>=1 else R(0)
older=coeff[k-2] if k>=2 else R(0)
factor=k*(k+1)+q-chi
divisor=2*(k+1)**2
coeff.append((factor*a-2*q*previous+q*older)/divisor)
if with_derivative:
da=derivative[k]
dp=derivative[k-1] if k>=1 else R(0)
do=derivative[k-2] if k>=2 else R(0)
derivative.append((factor*da-a-2*q*dp+q*do)/divisor)
return coeff,derivative
def boundary(chi,q,parity,terms):
coeff,derivative=endpoint_coefficients(chi,q,terms,True)
if parity:
return sum(coeff),sum(derivative)
return (sum(k*a for k,a in enumerate(coeff)),
sum(k*a for k,a in enumerate(derivative)))
def evaluate(coeff,x):
t=1-x
value=coeff[-1]
dt=x.parent()(0)
for a in reversed(coeff[:-1]):
dt=dt*t+value
value=value*t+a
return value,-dt
def finite(x):
return not x.is_infinity() and not x.is_NaN()
@lru_cache(maxsize=None)
def endpoint_function(S,n,dps,terms):
R=RealField(numberdb.bits(dps,losing=16))
q=(R.pi()*R(S)/2)**2
chi=R(str(pro_cv(0,n,float(R.pi()*S/2))))
for iteration in range(12):
residual,derivative=boundary(chi,q,n%2,terms)
assert finite(residual) and finite(derivative) and derivative
change=residual/derivative
chi-=change
if abs(change)<R(10)**(-(dps-30)):
break
else:
raise ArithmeticError('Independent eigenvalue iteration did not converge')
coeff,_=endpoint_coefficients(chi,q,terms)
return chi,tuple(coeff)
def independent_root(coeff,seed,dps):
R=coeff[0].parent()
x=R(str(seed))
for _ in range(15):
value,derivative=evaluate(coeff,x)
assert finite(value) and finite(derivative) and derivative
change=value/derivative
x-=change
if abs(change)<R(10)**(-(dps-35)):
assert 0<x<1
return x
raise ArithmeticError('Independent root iteration did not converge')
def agreement_digits(a,b):
R=RealField(1000)
error=abs(R(a)-R(b))/abs(R(b))
return 290 if error==0 else int(-error.log10())
def parse_decimal(text):
mantissa,marker,exponent=text.lower().partition('e')
exponent=int(exponent) if marker else 0
integer,dot,fraction=mantissa.partition('.')
assert dot and integer=='0'
unit=QQ(10)**(exponent-len(fraction))
return QQ(int(integer+fraction,10))*unit,unit
def controls():
# These nonzero values make a checker that returns zero everywhere fail.
c0,_=endpoint_coefficients(QQ(0),QQ(0),20)
c1,_=endpoint_coefficients(QQ(2),QQ(0),20)
c2,_=endpoint_coefficients(QQ(6),QQ(0),20)
assert evaluate(c0,QQ(1)/2)==(QQ(1),QQ(0))
assert evaluate(c1,QQ(1)/2)==(QQ(1)/2,QQ(1))
assert evaluate(c2,QQ(1)/2)==(-QQ(1)/8,QQ(3)/2)
R=RealField(900)
x=1/R(3).sqrt()
p2,_=endpoint_coefficients(R(6),R(0),20)
root=independent_root(p2,'0.58',220)
assert abs(root-x)<R(10)**(-180)
assert evaluate(p2,x-R('1e-100'))[0]*evaluate(p2,x+R('1e-100'))[0]<0
# At c=0 these are the exact Legendre polynomials, in both representations.
for n in range(12):
coeff,_=endpoint_coefficients(QQ(n*(n+1)),QQ(0),30)
assert all(a==0 for a in coeff[n+1:])
residual,derivative=boundary(QQ(n*(n+1)),QQ(0),n%2,30)
assert residual==0 and derivative!=0
leg=[QQ(0)]*max(2,n+1)
leg[n]=QQ(1)
for xq in (QQ(0),QQ(1)/3,-QQ(1)/2,QQ(1)):
assert evaluate(coeff,xq)==g.legendre_value(leg,xq)
return {'nonzero_endpoint_controls':3,'known_root':'P2 positive zero = 1/sqrt(3)',
'legendre_limit_modes':12,'exact_legendre_value_and_derivative_comparisons':48}
def stored_values():
tree=numberdb.table('T534')
result={}
nums=tree['Numbers']
if isinstance(nums,list):
for item in nums:
p=item['params']
result[tuple(int(p[k]) for k in ('m','S','n','j'))]=item['number']
else:
def walk(node,keys):
if len(keys)==4:
result[tuple(map(int,keys))]=node['number'] if isinstance(node,dict) else node
return
node=node.get('numbers',node)
for key,value in node.items():
walk(value,keys+[key])
walk(nums,[])
expected={tuple(int(p[k]) for k in g.ProlateZeros.parameters) for p in g.ProlateZeros().enumerate()}
assert set(result)==expected
return result
def run(use_stored=False):
start=time.monotonic()
base=Path(__file__).resolve().parent
report={'table':'T534','controls':controls(),'settings':g.SETTINGS,
'endpoint_settings':ENDPOINT_SETTINGS,'rows':[],
'read_from_database':use_stored,
'skill_manifest':'agents/runs/20261027T020220Z-build-skills/manifest.json',
'archived_skill_versions':json.loads((base/'skill-manifest.json').read_text()),
'source_sha256':{p:hashlib.sha256((base/p).read_bytes()).hexdigest()
for p in ('generate.py','spectral.py','checks.py','argument-selection.md')},
'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
R=RealField(1000)
for S in g.BANDWIDTHS:
for n in g.MODES:
roots=[g.roots(str(S),n,*setting) for setting in g.SETTINGS]
endpoints=[endpoint_function(S,n,*setting) for setting in ENDPOINT_SETTINGS]
# Normalize the endpoint solution against the independently normalized
# spectral expansion at x=1; compare actual nonzero values and derivatives.
spectral=g.eigenfunction(str(S),n,220,180)
at_one=g.legendre_value(spectral,spectral[0].parent()(1))[0]
for chi,coeff in endpoints:
F=chi.parent()
x=F(1)/3
value,derivative=evaluate(coeff,x)
reflected=evaluate(coeff,-x)[0]
assert abs(reflected-(-1)**n*value)<F('1e-150')*max(F(1),abs(value))
if n%2:
assert abs(evaluate(coeff,F(0))[0])<F('1e-150')
gv,gd=g.legendre_value(spectral,spectral[0].parent()(x))
assert abs(R(value)-R(gv/at_one))<R('1e-130')*max(R(1),abs(R(value)))
assert abs(R(derivative)-R(gd/at_one))<R('1e-130')*max(R(1),abs(R(derivative)))
prior=None
for j in range(1,n//2+1):
ref=R(roots[-1][j-1])
agreement=min(agreement_digits(rr[j-1],ref) for rr in roots[:-1])
assert agreement>=110
indep=[independent_root(coeff,g.coarse_roots(str(S),n)[j-1],setting[0])
for (_,coeff),setting in zip(endpoints,ENDPOINT_SETTINGS)]
independent_agreement=min(agreement_digits(ref,x) for x in indep)
assert independent_agreement>=110
p={'m':'0','S':str(S),'n':str(n),'j':str(j)}
text=stored[(0,S,n,j)] if stored is not None else g.ProlateZeros().value(p,100)
midpoint,ulp=parse_decimal(text)
assert 0<midpoint-ulp<midpoint+ulp<1
assert prior is None or prior<midpoint-ulp
prior=midpoint+ulp
assert abs(R(midpoint)-R(indep[-1]))<R(ulp)
sign_tests=[]
for (chi,coeff),setting in zip(endpoints,ENDPOINT_SETTINGS):
F=chi.parent()
lo=evaluate(coeff,F(midpoint-ulp))[0]
hi=evaluate(coeff,F(midpoint+ulp))[0]
assert finite(lo) and finite(hi) and lo*hi<0
sign_tests.append([int(lo.sign()),int(hi.sign())])
F=spectral[0].parent()
lo=g.legendre_value(spectral,F(midpoint-ulp))[0]
hi=g.legendre_value(spectral,F(midpoint+ulp))[0]
assert lo*hi<0
derivative=evaluate(endpoints[-1][1],indep[-1])[1]
assert derivative and finite(derivative)
# SciPy is an external low-precision check, not the source of the digits.
scipy_root=g.coarse_roots(str(S),n)[j-1]
scipy_error=abs(R(str(scipy_root))-ref)
assert scipy_error<R('1e-10')
report['rows'].append({'m':0,'S':S,'n':n,'j':j,'number':text,
'precision_truncation_agree_digits':agreement,
'independent_agree_digits':independent_agreement,
'endpoint_self_agree_digits':agreement_digits(*indep),
'endpoint_interval_signs':sign_tests,
'spectral_interval_signs':[int(lo.sign()),int(hi.sign())],
'endpoint_derivative':str(derivative.n(40)),
'scipy_root_absolute_error':str(scipy_error.n(40)),
'rounding_error_ulps':str((abs(R(midpoint)-R(indep[-1]))/R(ulp)).n(40))})
print('checked all positive zeros at S=%d'%S,file=sys.stderr,flush=True)
assert len(report['rows'])==120
report['minimum_precision_truncation_agreement']=min(r['precision_truncation_agree_digits'] for r in report['rows'])
report['minimum_independent_agreement']=min(r['independent_agree_digits'] for r in report['rows'])
report['minimum_endpoint_self_agreement']=min(r['endpoint_self_agree_digits'] for r in report['rows'])
report['maximum_scipy_absolute_error']=max(float(r['scipy_root_absolute_error']) for r in report['rows'])
report['root_count_indexing']='Every mode has floor(n/2) disjoint positive sign-change intervals. Parity, the exact odd central zero and the DLMF theorem of exactly n interior zeros identify all indices.'
report['limitations']='Finite-precision Taylor and Legendre truncation comparisons supply no certified infinite-problem error bound. All numerical signs are heuristic.'
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))