back to table · edit · history · where entries came from · files · download
5585 bytes, as of the version from 2026-10-27 02:15 (current). Recorded here, not run.
"""Zeros of prolate spheroidal wave functions -- numberdb.org/T534.
Run with SageMath and the NumberDB client; keep the attachments together:
sage -pip install numberdb
sage -python generate.py # verify every stored entry
sage -python generate.py --fill # fill an authorized draft
sage -python checks.py # independent numerical checks
sage -python checks.py --stored # check the stored decimal intervals
Filling sends numbers and attachments, without offering or publishing the
private draft. Authenticated operations use NUMBERDB_API_KEY, or a key piped
on stdin with NUMBERDB_KEY_FROM_STDIN=1. The campaign wrapper requests filling
with NUMBERDB_PUBLISH=1. The selected inputs are in argument-selection.md.
The generator uses normalized Legendre polynomials. Backwards recurrence
constructs the minimal eigenvector for each eigenvalue obtained by indexed
Sturm bisection. Newton iteration refines separately bracketed positive roots.
Arithmetic precision and Legendre truncation vary independently. At least 110
relative decimal digits must agree before storing 100. There is no certified
infinite-series remainder bound; the rigour is heuristic agreement-checked.
"""
import os
import sys
from functools import lru_cache
import numberdb.sage as numberdb
from sage.rings.rational_field import QQ
from scipy.special import pro_ang1
from scipy.optimize import brentq
from math import pi
import spectral as sp
BANDWIDTHS = (1, 2, 4, 8)
MODES = range(2, 12)
SETTINGS = ((160, 80), (220, 80), (220, 120), (220, 180))
REQUIRED_AGREEMENT = 110
def legendre_value(coefficients, x):
"""Evaluate sum a_k P_k(x) and its derivative by the three-term recurrence."""
R=x.parent()
p0,p1,d0,d1=R(1),x,R(0),R(1)
total=coefficients[0]+coefficients[1]*x
derivative=coefficients[1]
for k in range(1,len(coefficients)-1):
p2=((2*k+1)*x*p1-k*p0)/(k+1)
d2=((2*k+1)*(p1+x*d1)-k*d0)/(k+1)
total+=coefficients[k+1]*p2
derivative+=coefficients[k+1]*d2
p0,p1,d0,d1=p1,p2,d1,d2
return total,derivative
@lru_cache(maxsize=None)
def eigenfunction(S,n,dps,size):
R=sp.field(dps)
q=(R.pi()*R(QQ(S))/2)**2
chi=sp.spectrum(S,dps,size)[n]
diagonal,squares=sp.matrix(q,n%2,size)
off=[b.sqrt() for b in squares]
v=[R(0)]*size
v[-1]=R(1)
for i in range(size-2,-1,-1):
upper=off[i+1]*v[i+2] if i+2<size else R(0)
v[i]=((chi-diagonal[i+1])*v[i+1]-upper)/off[i]
norm=sum(b*b for b in v).sqrt()
v=[b/norm for b in v]
a=[R(0)]*(2*size)
for i,b in enumerate(v):
k=2*i+n%2
a[k]=b*(R(2*k+1)/2).sqrt()
if sum(a)<0:
a=[-b for b in a]
residual=abs((diagonal[0]-chi)*v[0]+off[0]*v[1])
assert residual<R(10)**(-(dps-20))
return tuple(a)
@lru_cache(maxsize=None)
def coarse_roots(S,n):
"""Double precision supplies seeds only; the expected count is compulsory."""
c=pi*float(QQ(S))/2
f=lambda x:pro_ang1(0,n,c,x)[0]
roots=[]
a=1e-8
fa=f(a)
for k in range(1,1001):
b=k/1001
fb=f(b)
if fa*fb<0:
roots.append(brentq(f,a,b,xtol=1e-15))
a,fa=b,fb
assert len(roots)==n//2, (S,n,roots)
return tuple(roots)
@lru_cache(maxsize=None)
def roots(S,n,dps,size):
R=sp.field(dps)
a=eigenfunction(S,n,dps,size)
result=[]
for seed in coarse_roots(S,n):
x=R(str(seed))
for iteration in range(15):
value,derivative=legendre_value(a,x)
assert derivative
change=value/derivative
x-=change
if abs(change)<R(10)**(-(dps-20)):
break
else:
raise ArithmeticError('Root iteration did not converge')
assert 0<x<1
assert not result or x>result[-1]
result.append(x)
assert len(result)==n//2
return tuple(result)
class ProlateZeros(numberdb.Generator):
table='T534'
parameters=('m','S','n','j')
type='R'
digits=100
rigour='heuristic (agreement-checked)'
files=('generate.py','spectral.py','checks.py','argument-selection.md',
'evidence.json','skill-manifest.json')
def enumerate(self):
for S in BANDWIDTHS:
for n in MODES:
for j in range(1,n//2+1):
yield {'m':'0','S':str(QQ(S)),'n':str(n),'j':str(j)}
def value(self,params,digits):
if int(params['m'])!=0 or digits!=100:
raise ValueError('Only order zero and measured 100-digit settings are supported')
S,n,j=str(params['S']),int(params['n']),int(params['j'])
if n not in MODES or not 1<=j<=n//2:
raise ValueError('Unsupported mode or zero index')
values=[roots(S,n,*setting)[j-1] for setting in SETTINGS]
R=sp.field(220)
ref=R(values[-1])
assert all(abs(R(v)-ref)<abs(ref)*R(10)**(-REQUIRED_AGREEMENT) for v in values[:-1])
return sp.decimal_string(ref,digits)
if __name__=='__main__':
if os.environ.get('NUMBERDB_KEY_FROM_STDIN')=='1':
os.environ['NUMBERDB_API_KEY']=sys.stdin.read().strip()
generator=ProlateZeros()
if '--fill' in sys.argv or os.environ.get('NUMBERDB_PUBLISH')=='1':
for p in generator.enumerate():
generator.value(p,generator.digits)
print(generator.publish(message='120 positive angular zeros checked against an independent endpoint ODE'))
else:
report=generator.verify(sample=None)
print(report)
sys.exit(0 if report.ok else 1)