back to table · edit · history · where entries came from · files · download
8828 bytes, as of the version from 2026-10-27 09:57 (current). Recorded here, not run.
"""Reproduce checks for Lamé band-edge eigenvalues -- numberdb.org/T544.
Download generate.py, check.py, and skill-manifest.json beside this file.
Run: sage -python checks.py > numerical-evidence.json
Set CHECK_STORED=1 and NUMBERDB_API_KEY for the draft readback check.
If a key is piped on stdin, set NUMBERDB_KEY_FROM_STDIN=1 instead.
Exact certification uses rational Sturm counts. Independent numerical checking
uses Fourier Galerkin matrices of the original periodic differential operator,
with NumPy/SciPy and no sector-polynomial code. This comparison is a double-
precision diagnostic, not the basis of the 100-digit enclosure claim.
"""
import os,sys,json,hashlib
from pathlib import Path
from fractions import Fraction
from decimal import Decimal
import numberdb.sage as numberdb
from sage.rings.rational_field import QQ
from sage.rings.real_mpfi import RealIntervalField
from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
from sage.version import version as sage_version
import numpy as np
import scipy
from scipy.special import ellipj,ellipk
from scipy.linalg import eigvalsh
from numberdb._write import to_text
from numberdb._generate import _stored_entries
from generate import (LameBandEdges,band_edges,spectral_polynomial,root_certificates,
sturm_sequence,variations,sector_polynomials,S,h)
from check import exactness,measure,agrees_with,names_its_rings
HERE=Path(__file__).resolve().parent
def text_interval(text):
text=str(text)
q=Fraction(text)
center=QQ(q.numerator)/q.denominator
if '.' not in text and 'e' not in text.lower():
return center,center
exp=Decimal(text).as_tuple().exponent
ulp=QQ(10)**exp
return center-ulp,center+ulp
def operator_edges(n,m,modes):
"""Fourier Galerkin on [0,2K], both periodic and antiperiodic sectors."""
m=float(m)
K=ellipk(m)
N=2048
x=2*K*np.arange(N)/N
V=n*(n+1)*m*ellipj(x,m)[0]**2
coeff=np.fft.fft(V)/N
modes_array=np.arange(-modes,modes+1)
potential=coeff[(modes_array[:,None]-modes_array[None,:]) % N]
assert np.max(np.abs(potential-potential.conj().T))<1e-11
eigenvalues=[]
for shift in (0.0,0.5):
kinetic=((modes_array+shift)*np.pi/K)**2
vals=eigvalsh(potential+np.diag(kinetic))
eigenvalues.extend(vals[:2*n+2])
return sorted(eigenvalues)[:2*n+1]
def maier_polynomial(n,m):
"""Independent transcription of Maier equations (54)-(56)."""
if n==1:
return (h-1)*(h-m)*(h-m-1)
if n==2:
return (h*h-4*(m+1)*h+12*m)*(h-m-1)*(h-4*m-1)*(h-m-4)
if n==3:
return ((h-4*m-4)*(h*h-2*(2*m+5)*h+3*(8*m+3))
*(h*h-2*(5*m+2)*h+3*(3*m*m+8*m))
*(h*h-10*(m+1)*h+3*(3*m*m+26*m+3)))
raise ValueError(n)
def main():
if os.environ.get('NUMBERDB_KEY_FROM_STDIN')=='1':
key=sys.stdin.read().strip()
if '=' in key and key.split('=',1)[0].isidentifier():
key=key.split('=',1)[1].strip().strip('\"\'')
numberdb.configure(api_key=key)
gen=LameBandEdges()
# Run numerical controls BEFORE accepting the operator comparison.
control_errors=[]
for n in (1,6,12):
got=operator_edges(n,0,48)
known=[0]+[k*k for k in range(1,n+1) for _ in range(2)]
err=max(abs(a-b) for a,b in zip(got,known))
assert err<1e-11
control_errors.append(err)
for m in (QQ(1)/4,QQ(1)/2,QQ(3)/4):
got=operator_edges(1,m,48)
known=[float(m),1.0,float(1+m)]
err=max(abs(a-b) for a,b in zip(got,known))
assert err<1e-10
control_errors.append(err)
assert names_its_rings(str(HERE/'generate.py'))==[]
params=list(gen.enumerate())
assert len(params)==264
raw={tuple(p.values()):gen.value(p,100) for p in params}
assert exactness(raw)==[]
rows=[{'params':p,'number':to_text(raw[tuple(p.values())],100)} for p in params]
written={tuple(r['params'].values()):r['number'] for r in rows}
mode='local'
if os.environ.get('CHECK_STORED')=='1':
doc=numberdb.table('T544')
saved=_stored_entries(doc)
assert len(saved)==264
written={k:saved[','.join(map(str,k))] for k in raw}
assert list(saved)==[','.join(map(str,p.values())) for p in params]
mode='database readback'
checks=[]
certificates=[]
max_error=0.0
rational_count=0
widest=QQ(0)
min_gap=None
maier_cases={}
groups=[]
for p in params:
nm=(p['n'],p['m'])
if nm not in groups: groups.append(nm)
# Exact matches to external low-order formulas.
for n,mtext in groups:
m=QQ(mtext)
if n<=3: maier_cases[(n,mtext)]=spectral_polynomial(n,mtext)
assert agrees_with(maier_cases,lambda key:maier_polynomial(key[0],QQ(key[1])))==[]
corpus=numberdb.table('T543')
BM=PolynomialRing(QQ,('h','m'))
# T543 is an additional corpus consistency check, not independence evidence.
corpus_checks=0
for n,mtext in groups:
m=QQ(mtext)
poly=spectral_polynomial(n,mtext)
seq=sturm_sequence(poly)
assert poly == -spectral_polynomial(n,str(1-m))(n*(n+1)-h)
if str(n) in corpus['Numbers']:
value=corpus['Numbers'][str(n)]
if isinstance(value,dict): value=value['number']
assert poly == S(BM(value)(h,m))
corpus_checks+=1
if m==QQ(1)/2:
assert poly(QQ(n*(n+1))/2)==0
assert band_edges(n,mtext)[n]==QQ(n*(n+1))/2
direct48=operator_edges(n,m,48)
direct64=operator_edges(n,m,64)
intervals=[]
certs=root_certificates(n,mtext)
errors=[]
for j,(lo,hi,f) in enumerate(certs):
text=written[(n,mtext,j)]
wl,wu=text_interval(text)
assert wl<=lo<=hi<=wu
if wl==wu:
assert poly(wl)==0 and f.degree()==1 and wl==lo==hi
rational_count+=1
else:
assert variations(seq,wl)-variations(seq,wu)==1
assert poly(wl)*poly(wu)<0
widest=max(widest,hi-lo)
if intervals:
gap=wl-intervals[-1][1]
assert gap>0
min_gap=gap if min_gap is None else min(min_gap,gap)
intervals.append((wl,wu))
midpoint=float((lo+hi)/2)
errors.append(max(abs(direct48[j]-midpoint),abs(direct64[j]-midpoint)))
assert errors[-1]<2e-10
# Complementary parameter relation on written values, all rows.
cl,cu=text_interval(written[(n,str(1-m),2*n-j)])
assert wl+cl<=n*(n+1)<=wu+cu
certificates.append({'params':{'n':n,'m':mtext,'j':j},'number':text,
'factor':str(f),'lo':str(lo),'hi':str(hi)})
max_error=max(max_error,max(errors))
checks.append({'n':n,'m':mtext,'edge_count':2*n+1,
'max_operator_absolute_error':max(errors),
'max_48_64_mode_difference':max(abs(a-b) for a,b in zip(direct48,direct64))})
result={'table':'T544','mode':mode,'skill_manifest':'agents/runs/20261027T094619Z-build-skills/manifest.json',
'skill_versions':json.loads((HERE/'skill-manifest.json').read_text()),
'versions':{'sage':sage_version,'numberdb':numberdb.__version__,'numpy':np.__version__,'scipy':scipy.__version__},
'source_sha256':{name:hashlib.sha256((HERE/name).read_bytes()).hexdigest() for name in ['generate.py','checks.py','check.py']},
'entries':len(rows),'exact_rational_entries':rational_count,'irrational_entries':len(rows)-rational_count,
'rendered_measure':measure(written),'enclosure_width_max':str(widest),
'enclosure_width_max_float':float(widest),'minimum_written_interval_gap':float(min_gap),
'controls':{'known_free_operator_orders':[1,6,12],'known_n1_parameters':['1/4','1/2','3/4'],
'maximum_error':max(control_errors)},
'exact_checks':{'Maier_54_56_cases':len(maier_cases),'T543_specializations':corpus_checks,
'complementary_polynomial_cases':len(groups),'lemniscatic_central_roots':12,
'written_interval_root_and_order_certificates':len(rows),
'written_complementary_identity_checks':len(rows)},
'independent_operator_check':{'method':'Fourier Galerkin of original operator, 48 and 64 modes per side, 2048 sampled potential points, both boundary multipliers',
'certifies_100_digits':False,'max_absolute_error':max_error,'cases':checks},
'certificates':certificates,'all_checks_passed':True}
print(json.dumps(result,indent=2))
if __name__=='__main__':main()