back to table · edit · history · where entries came from · files · download
6919 bytes, as of the version from 2026-10-26 16:08. Recorded here, not run.
"""Independent exact checks for numberdb.org/T520.
sage -python checks.py # check the generator locally
sage -python checks.py --stored # check every stored coefficient
Download generate.py beside this file. The stored check reads the draft with
NUMBERDB_API_KEY, or reads a key from stdin with NUMBERDB_KEY_FROM_STDIN=1.
No values, documents, review states or publication states are written.
This file implements the differential equation in monomials, without calling
the generator's perturbation or Legendre-multiplication code to form expected
values. DLMF 30.3.9--12 provide external formulas for k=0,...,4.
"""
import numberdb.sage as numberdb
from sage.rings.rational_field import QQ
from functools import lru_cache
from math import factorial
import json
import os
import sys
from generate import SpheroidalCoefficients, N_MAX, K_MAX, perturbation
def rat(a, b=1):
return QQ(a) / QQ(b)
@lru_cache(maxsize=None)
def monomial_series(n, order):
"""Solve in powers of x, with coefficient of x^n fixed to one."""
d = n*(n+1)
g0 = {n: QQ(1)}
for l in range(n-2, -1, -2):
g0[l] = rat((l+2)*(l+1), l*(l+1)-d)*g0[l+2]
gs, aa = [g0], [QQ(d)]
for r in range(1, order+1):
rhs = {l+2: -v for l,v in gs[-1].items()}
for j in range(1,r):
for l,v in gs[r-j].items():
rhs[l] = rhs.get(l,QQ(0)) + aa[j]*v
current = {}
for l in range(n+2*r,n,-2):
current[l] = (rhs.get(l,QQ(0))+(l+2)*(l+1)*current.get(l+2,QQ(0))) / QQ(l*(l+1)-d)
ar = -(n+2)*(n+1)*current.get(n+2,QQ(0))-rhs.get(n,QQ(0))
aa.append(ar)
current[n] = QQ(0)
for l in range(n-2,-1,-2):
current[l] = (rhs.get(l,QQ(0))+ar*g0[l]+(l+2)*(l+1)*current[l+2]) / QQ(l*(l+1)-d)
gs.append(current)
return tuple(aa)
def dlmf(n):
"""DLMF 30.3.9--12 at m=0, shifting ell_2 by one."""
l0 = QQ(n*(n+1))
l2 = (-1 + rat(1,(2*n-1)*(2*n+3)))/QQ(2)
lower=(n-1)**2*n**2
upper=(n+1)**2*(n+2)**2
l4=(rat(lower,(2*n-3)*(2*n-1)**3*(2*n+1))
-rat(upper,(2*n+1)*(2*n+3)**3*(2*n+5)))/QQ(2)
l6=-(rat(upper,(2*n-1)*(2*n+1)*(2*n+3)**5*(2*n+5)*(2*n+7))
-rat(lower,(2*n-5)*(2*n-3)*(2*n-1)**5*(2*n+1)*(2*n+3)))
A=(rat(lower,(2*n-5)**2*(2*n-3)*(2*n-1)**7*(2*n+1)*(2*n+3)**2)
-rat(upper,(2*n-1)**2*(2*n+1)*(2*n+3)**7*(2*n+5)*(2*n+7)**2))
B=(rat(((n-3)*(n-2)*(n-1)*n)**2,(2*n-7)*(2*n-5)**2*(2*n-3)**3*(2*n-1)**4*(2*n+1))
-rat(((n+1)*(n+2)*(n+3)*(n+4))**2,(2*n+1)*(2*n+3)**4*(2*n+5)**3*(2*n+7)**2*(2*n+9)))
C=(rat(upper**2,(2*n+1)**2*(2*n+3)**7*(2*n+5)**2)
-rat(lower**2,(2*n-3)**2*(2*n-1)**7*(2*n+1)**2))
D=rat(lower*upper,(2*n-3)*(2*n-1)**4*(2*n+1)**2*(2*n+3)**4*(2*n+5))
l8=2*A+B/QQ(16)+C/QQ(8)+D/QQ(2)
return (l0, l2+1, l4, l6, l8)
@lru_cache(maxsize=None)
def legendre_rodrigues(l):
return {l-2*j: rat((-1)**j * factorial(2*l-2*j),
2**l * factorial(j)*factorial(l-j)*factorial(l-2*j))
for j in range(l//2+1)}
def add_scaled(result, poly, scale):
for j,c in poly.items():
result[j]=result.get(j,QQ(0))+scale*c
def residual_check(n, aa, vectors):
polynomials=[]
for vector in vectors:
poly={}
for l,c in vector.items():
add_scaled(poly,legendre_rodrigues(l),c)
polynomials.append(poly)
for r, poly in enumerate(polynomials):
residual={}
for l,c in poly.items():
residual[l]=residual.get(l,QQ(0))+(l*(l+1)-aa[0])*c
if l>=2: residual[l-2]=residual.get(l-2,QQ(0))-l*(l-1)*c
if r:
add_scaled(residual,{l+2:c for l,c in polynomials[r-1].items()},QQ(1))
for j in range(1,r+1):
add_scaled(residual,polynomials[r-j],-aa[j])
assert all(c==0 for c in residual.values()), ('ODE residual',n,r)
def controls():
# These controls must return nonzero known answers, not just zero residuals.
assert legendre_rodrigues(0)=={0:QQ(1)}
assert legendre_rodrigues(1)=={1:QQ(1)}
assert legendre_rodrigues(2)=={2:rat(3,2),0:rat(-1,2)}
assert monomial_series(0,2)==(QQ(0),rat(1,3),rat(-2,135))
assert monomial_series(1,1)==(QQ(2),rat(3,5))
assert dlmf(0)==(QQ(0),rat(1,3),rat(-2,135),rat(4,8505),rat(26,1913625))
def check_values(values):
controls()
expected={(0,n,k) for n in range(N_MAX+1) for k in range(K_MAX+1)}
assert set(values)==expected, 'Enumeration differs from the stated selection'
dlmf_count=0
for n in range(N_MAX+1):
aa=tuple(values[(0,n,k)] for k in range(K_MAX+1))
assert aa==monomial_series(n,K_MAX), ('monomial disagreement',n)
for k,external in enumerate(dlmf(n)):
assert aa[k]==external, ('DLMF disagreement',n,k)
dlmf_count+=1
# Use the supplied values, not recomputed eigenvalues, in this residual.
residual_check(n,aa,perturbation(n,K_MAX)[1])
return {'entries':len(values),'monomial_matches':len(values),
'DLMF_formula_matches':dlmf_count,'exact_ODE_residuals':len(values),
'controls':'passed','selection':'m=0, 0<=n<=11, 0<=k<=12'}
def parse_rational(value):
# Avoid QQ(text): modular Sage may reach for uninitialised real_mpfr.
parts = str(value).split('/')
if len(parts) == 1:
return QQ(int(parts[0]))
if len(parts) == 2:
return QQ(int(parts[0])) / QQ(int(parts[1]))
raise ValueError('Expected an exact integer or rational')
def flatten(document):
"""Read either NumberDB nested Numbers or flat records; preserve all keys."""
numbers=document.get('Numbers',[])
if isinstance(numbers,list):
for row in numbers:
yield tuple(int(row['params'][p]) for p in ('m','n','k')),parse_rational(row['number'])
else:
def walk(block,depth,keys):
if depth==3:
yield tuple(keys),parse_rational(block['number'] if isinstance(block,dict) else block)
return
for key,child in block.items():
if key in ('comment','param-latex'): continue
yield from walk(child,depth+1,keys+[int(key)])
yield from walk(numbers,0,[])
if __name__=='__main__':
if os.environ.get('NUMBERDB_KEY_FROM_STDIN')=='1':
secret=sys.stdin.read().strip()
if '=' in secret and secret.split('=',1)[0].isidentifier():
secret=secret.split('=',1)[1].strip().strip('\"').strip("'")
os.environ['NUMBERDB_API_KEY']=secret
if '--stored' in sys.argv or os.environ.get('T520_CHECK_STORED')=='1':
values=dict(flatten(numberdb.table('T520')))
source='live database'
else:
g=SpheroidalCoefficients()
values={tuple(p[x] for x in g.parameters):g.value(p,100) for p in g.enumerate()}
source='local generator'
print(json.dumps({'source':source,**check_values(values)},indent=2))