"""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()