back to table · edit · history · where entries came from · files · download
9046 bytes, as of the version from 2026-10-30 13:35 (current). Recorded here, not run.
"""Independent numerical and graph checks for numberdb.org/T626.
Run: sage -python checks.py > numerical-evidence.json
Dependencies: numberdb, mpmath, numpy, and the attached generate.py.
Quadrature/finite tori are independent checks, not the source of certification.
For a database readback, run stored_checks.py after filling the draft.
"""
import json
import sys
from collections import Counter, deque
from decimal import Decimal
from functools import lru_cache
from fractions import Fraction
from pathlib import Path
import numberdb.sage as numberdb
from numberdb._write import to_text
from sage.rings.rational_field import QQ
import mpmath as mp
import numpy as np
import generate as g
def key(p):
return ','.join(str(p[k]) for k in g.Resistance.parameters)
def unbox(v):
return v['number'] if isinstance(v,dict) else v
def endpoint_rationals(v):
if hasattr(v,'lower'): return v.lower().exact_rational(),v.upper().exact_rational()
return QQ(v),QQ(v)
def written_bounds(text):
if '.' not in text and 'e' not in text.lower():
q=QQ(text)
return q,q
d=Decimal(text)
center=QQ(Fraction(d))
unit=QQ(10)**d.as_tuple().exponent
return center-unit,center+unit
def bfs(lattice,radius):
start=('A',0,0)
dist={start:0}
queue=deque([start])
while queue:
s,m,n=queue.popleft()
if dist[s,m,n] == radius: continue
if lattice == 'triangular':
neighbors=[('A',m+a,n+b) for a,b in [(1,0),(-1,0),(0,1),(0,-1),(1,1),(-1,-1)]]
elif s == 'A': neighbors=[('B',m,n),('B',m-1,n),('B',m-1,n-1)]
else: neighbors=[('A',m,n),('A',m+1,n),('A',m+1,n+1)]
for v in neighbors:
if v not in dist:
dist[v]=dist[s,m,n]+1
queue.append(v)
return dist
@lru_cache(None)
def integral(m,n,dps):
with mp.workdps(dps):
diff=abs(m-n)
total=m+n
def f(y):
if not y: return mp.mpf(diff)
c=mp.cos(y)
if y == mp.pi/2:
return (1-mp.cos(total*y))/2 if diff == 0 else mp.mpf(1)/2
x=2*mp.asinh(mp.sin(y/2)*mp.sqrt((c+2)/c))
decay=mp.exp(-diff*x)
numerator=-mp.expm1(-diff*x)+2*decay*mp.sin(total*y/2)**2
return numerator/(mp.sin(y)*mp.sqrt(3+mp.sin(y)**2))
return mp.quad(f,[0,mp.pi/8,mp.pi/4,3*mp.pi/8,mp.pi/2],method='gauss-legendre')/mp.pi
def quadrature(p,dps):
m,n=int(p['m']),int(p['n'])
with mp.workdps(dps):
if p['lattice']=='triangular': return integral(m,n,dps)
if p['sublattice']=='A': return 3*integral(m,n,dps)
return integral(m,n,dps)+integral(m+1,n,dps)+integral(m+1,n+1,dps)
def torus(N):
theta=2*np.pi*np.arange(N)/N
x,y=theta[:,None],theta[None,:]
# Build each network's Fourier Laplacian from its own declared edges.
lam=6-2*np.cos(x)-2*np.cos(y)-2*np.cos(x+y)
lam[0,0]=1
inv=1/lam
inv[0,0]=0
green=np.fft.ifft2(inv).real
t=2*(green[0,0]-green)
coupling=1+np.exp(1j*x)+np.exp(1j*(x+y))
det=9-abs(coupling)**2
det[0,0]=1
diag=3/det
cross=coupling/det
# At zero wave vector the nonconstant basis mode has eigenvalue six.
diag[0,0]=1/12
cross[0,0]=-1/12
aa=np.fft.ifft2(diag).real
ba=np.fft.ifft2(cross).real
return t,2*(aa[0,0]-aa),2*(aa[0,0]-ba)
def main():
gen=g.Resistance()
params=list(gen.enumerate())
assert len(params)==len({key(p) for p in params})==84
counts=Counter((p['lattice'],p['sublattice']) for p in params)
assert counts=={('triangular','A'):48,('honeycomb','A'):15,('honeycomb','B'):21}
graph_summary={}
for lattice in ['triangular','honeycomb']:
dist=bfs(lattice,12)
reps=set()
for (s,m,n),d in dist.items():
if not d: continue
expected=g.distance(m,n) if lattice=='triangular' else g.honey_distance(s,m,n)
assert d==expected
representative=g.canonical(m,n) if s=='A' else max(g.b_orbit(m,n))
reps.add((s,*representative))
selected={(p['sublattice'],int(p['m']),int(p['n'])) for p in params if p['lattice']==lattice}
assert reps==selected
shells=dict(sorted(Counter(dist.values()).items()))
assert all(shells[d]==(6 if lattice=='triangular' else 3)*d for d in range(1,13))
graph_summary[lattice]={'vertices_including_origin':len(dist),'shell_counts':shells,'orbits':len(reps)}
# Published small values, evaluated before trusting quadrature.
mp.mp.dps=160
controls={(0,0):mp.mpf(0),(1,0):mp.mpf(1)/3,(2,0):mp.mpf(8)/3-4*mp.sqrt(3)/mp.pi,
(2,1):-mp.mpf(2)/3+2*mp.sqrt(3)/mp.pi,(3,1):-5+10*mp.sqrt(3)/mp.pi}
for mn,expected in controls.items():
assert abs(integral(*mn,125)-expected)<mp.mpf('1e-118'),(mn,'quadrature control')
print('Quadrature controls passed.',file=sys.stderr,flush=True)
records=[]
min_bits=10000
max_error=mp.mpf(0)
max_precision_difference=mp.mpf(0)
max_cancellation=mp.mpf(0)
for i,p in enumerate(params):
poly=g.coefficients(p)
assert all(c.parent() is QQ for c in poly.list())
v=unbox(gen.value(p,100))
lo,hi=endpoint_rationals(v)
text=to_text(v,100)
wlo,whi=written_bounds(text)
assert wlo<=lo<=hi<=whi,(key(p),'serialization')
if hi!=lo:
assert v.is_finite() and lo>0
relative=(hi-lo)/lo
accuracy=relative.denominator().nbits()-relative.numerator().nbits()-1
assert relative < QQ(2)**(-accuracy)
assert accuracy>=350,(key(p),accuracy)
min_bits=min(min_bits,int(accuracy))
else: accuracy=None
# Independent integral has no access to the recurrence coefficients.
q1,q2=quadrature(p,125),quadrature(p,150)
a,b=mp.mpf(str(poly[0])),mp.mpf(str(poly[1]))
expected=a+b*mp.sqrt(3)/mp.pi
error=abs(q2-expected)
difference=abs(q1-q2)
assert error<mp.mpf('1e-140'),(key(p),str(error))
assert difference<mp.mpf('1e-118'),(key(p),str(difference))
max_error=max(max_error,error)
max_precision_difference=max(max_precision_difference,difference)
cancellation=mp.log((abs(a)+abs(b*mp.sqrt(3)/mp.pi))/abs(expected),2)
max_cancellation=max(max_cancellation,cancellation)
records.append({'params':p,'number':text,'exact_coefficients':[str(poly[0]),str(poly[1])],
'relative_accuracy_bits_lower_bound':accuracy,'integral_150_digits':mp.nstr(q2,145),
'quadrature_absolute_error':mp.nstr(error,8)})
if (i+1)%12==0: print('Checked %d/84 entries.'%(i+1),file=sys.stderr,flush=True)
# Harmonicity at every selected triangular displacement, including the
# additional outer neighbors needed to check the boundary shell.
for p in params:
m,n=int(p['m']),int(p['n'])
if p['lattice']=='triangular':
assert 6*g.triangular(m,n)==sum(g.triangular(m+a,n+b) for a,b in g.STEPS)
else:
s=p['sublattice']
if s=='A': neighbors=[('B',m,n),('B',m-1,n),('B',m-1,n-1)]
else: neighbors=[('A',m,n),('A',m+1,n),('A',m+1,n+1)]
assert 3*g.coefficients(p)==sum(g.coefficients(dict(lattice='honeycomb',sublattice=s1,m=str(a),n=str(b))) for s1,a,b in neighbors)
torus_summary=[]
previous=None
for N in [64,128,256]:
t,ha,hb=torus(N)
assert abs(t[1,0]-(N*N-1)/(3*N*N))<1e-12
assert abs(hb[0,0]-(2*N*N-1)/(3*N*N))<1e-12
errors=[]
for p in params:
array=t if p['lattice']=='triangular' else ha if p['sublattice']=='A' else hb
observed=array[int(p['m'])%N,int(p['n'])%N]
target=float(unbox(gen.value(p,100)))
errors.append(abs(observed-target))
if previous is not None: assert all(b<a for a,b in zip(previous,errors))
previous=errors
torus_summary.append({'side':N,'maximum_absolute_error':max(errors),'adjacent_controls_passed':True})
evidence={'table':'T626','skill_manifest':'agents/runs/20261030T125612Z-build-skills/manifest.json',
'skill_sha256': '28ecbeeccb0c5e5d6eb541b5a5b0f88aec47b154374d909085116d34f5629199',
'entries':84,'group_counts':{','.join(k):v for k,v in counts.items()},'graph_checks':graph_summary,
'working_precision_bits':numberdb.bits(100,losing=g.WORKING_GUARD),'minimum_relative_accuracy_bits':min_bits,
'maximum_cancellation_bits':mp.nstr(max_cancellation,12),'quadrature_working_digits':[125,150],
'maximum_quadrature_error':mp.nstr(max_error,12),'maximum_precision_difference':mp.nstr(max_precision_difference,12),
'exact_harmonicity_checks':84,'serialized_enclosures_checked':84,'torus_checks':torus_summary,
'versions':{'numberdb':getattr(numberdb,'__version__','unknown'),'mpmath':mp.__version__,'numpy':np.__version__},
'records':records}
print(json.dumps(evidence,indent=2))
if __name__=='__main__': main()