back to table · edit · history · where entries came from · files · download
9379 bytes, as of the version from 2026-10-29 01:41 (current). Recorded here, not run.
"""Independent numerical and structural checks for numberdb.org/T596.
Run: sage -python checks.py [--stored stored-document.json]
The generator uses exact recurrences. This checker separately integrates
Cserti (27), twice, with mpmath tanh-sinh quadrature. Quadrature is an
agreement check, not the source of the certified error claims. It also
checks finite periodic networks by a discrete Fourier solve and a known
finite-torus adjacent-vertex control. All outputs are JSON.
"""
import json
import sys
import platform
import hashlib
from pathlib import Path
from decimal import Decimal
from fractions import Fraction
import numberdb.sage as numberdb
import numberdb as client_module
from numberdb._write import to_text
from sage.rings.rational_field import QQ
from sage.rings.real_mpfi import RealIntervalField
import sage.version
import mpmath as mp
import numpy as np
from generate import SquareLatticeResistance, coefficient_table, coefficients, certified_value
from check import exactness, measure, names_its_rings, agrees_with
HERE = Path(__file__).resolve().parent
def orbit(u, v):
return max(abs(u), abs(v)), min(abs(u), abs(v))
def mp_rational(q):
return mp.mpf(str(q.numerator())) / mp.mpf(str(q.denominator()))
def quadrature(m, n):
def integrand(y):
if not y:
return mp.mpf(m)
# cosh(s)=2-cos(y), with cancellation avoided at y=0.
s = 2*mp.asinh(mp.sin(y/2))
decay = mp.exp(-m*s)
numerator = -mp.expm1(-m*s) + 2*decay*mp.sin(n*y/2)**2
return numerator/mp.sinh(s)
return mp.quad(integrand, [0, mp.pi/4, mp.pi/2, mp.pi])/mp.pi
def text_interval(text):
if '/' in text or ('.' not in text and 'e' not in text.lower()):
q=Fraction(text)
return QQ(q.numerator)/QQ(q.denominator), QQ(q.numerator)/QQ(q.denominator)
q=Fraction(text)
mid=QQ(q.numerator)/QQ(q.denominator)
unit=QQ(10)**int(Decimal(text).as_tuple().exponent)
return mid-unit, mid+unit
def run():
gen=SquareLatticeResistance()
selected=[(int(p['m']),int(p['n'])) for p in gen.enumerate()]
# Derive the graph ball in signed lattice coordinates, without using enumerate.
expected=sorted({orbit(u,v) for u in range(-20,21) for v in range(-20,21)
if 0<abs(u)+abs(v)<=20},key=lambda q:(sum(q),q[1]))
assert selected==expected and len(selected)==120
counts={'axial':sum(n==0 for m,n in selected),
'diagonal':sum(m==n for m,n in selected),
'asymmetric':sum(0<n<m for m,n in selected)}
assert counts=={'axial':20,'diagonal':10,'asymmetric':90}
polys=coefficient_table(21)
assert not exactness(polys)
# Diagonal sum, derived by integrating along a diagonal, not its recurrence.
diagonal={k:polys[k,k] for k in range(1,11)}
variable=next(iter(diagonal.values())).parent().gen()
assert not agrees_with(diagonal,lambda k:2*variable*sum((QQ(1)/QQ(2*j-1) for j in range(1,k+1)),QQ(0)))
exact_harmonic=0
for m,n in selected:
neighbors=[orbit(m+1,n),orbit(m-1,n),orbit(m,n+1),orbit(m,n-1)]
assert 4*polys[m,n]==sum((polys[k] for k in neighbors),polys[0,0])
exact_harmonic+=1
values={k:certified_value(*k) for k in selected}
assert not exactness(values)
assert not names_its_rings(str(HERE/'generate.py'))
written={k:to_text(v,digits=100) for k,v in values.items()}
rows=[]; min_accuracy_bits=1000000; max_precision=0
for k,v in values.items():
low,high=text_interval(written[k])
if k==(1,0):
assert v==QQ(1)/2 and written[k]=='1/2'
else:
interval=RealIntervalField(v.parent().precision())(v)
l,h=(q.exact_rational() for q in interval.endpoints())
assert low<=l<=h<=high and l>0
ratio=l/(h-l)
conservative_bits=ratio.numerator().nbits()-ratio.denominator().nbits()-1
min_accuracy_bits=min(min_accuracy_bits,int(conservative_bits))
max_precision=max(max_precision,int(v.parent().precision()))
# Known answers run before using the integrator on the selected family.
controls=[]
independent={}
for dps in (120,145):
with mp.workdps(dps):
for pair,known in [((0,0),mp.mpf(0)),((1,0),mp.mpf(1)/2),((1,1),2/mp.pi)]:
computed=quadrature(*pair)
err=abs(computed-known)
assert err < mp.mpf(10)**(-dps+8)
controls.append({'digits':dps,'pair':pair,'absolute_error':mp.nstr(err,8)})
for k in selected:
q=quadrature(*k)
a,b=coefficients(*k)
exact_eval=mp_rational(a)+mp_rational(b)/mp.pi
assert abs(q-exact_eval)<mp.mpf(10)**(-dps+18)
independent[k,dps]=mp.nstr(q,dps)
max_residual=mp.mpf(0)
with mp.workdps(150):
for k in selected:
q=mp.mpf(independent[k,145]); earlier=mp.mpf(independent[k,120])
assert abs(q-earlier)<mp.mpf('1e-112')
low,high=text_interval(written[k])
if low==high:
assert abs(q-mp_rational(low))<mp.mpf('1e-135')
else:
assert mp_rational(low)<q<mp_rational(high)
a,b=coefficients(*k)
residual=abs(q-(mp_rational(a)+mp_rational(b)/mp.pi))
max_residual=max(max_residual,residual)
rows.append({'m':k[0],'n':k[1],'a':str(a),'b':str(b),'number':written[k],
'quadrature_120':independent[k,120],'quadrature_145':independent[k,145]})
# A separate finite-network solve: Fourier diagonalization of L=4I-A.
torus={}; torus_summary=[]
for size in (128,256,512):
c=np.cos(2*np.pi*np.arange(size)/size)
denominator=2-c[:,None]-c[None,:]
denominator[0,0]=1
kernel=1/denominator;kernel[0,0]=0
green=np.fft.ifft2(kernel).real
resistances=green[0,0]-green
control=abs(float(resistances[1,0])-(1-1/size**2)/2)
assert control<1e-13
errors=[]
for k in selected:
target=float(values[k])
error=target-float(resistances[k])
assert error>0
torus[size,k]=error;errors.append(error)
torus_summary.append({'side':size,'adjacent_control_error':control,'max_error':max(errors)})
ratios=[torus[256,k]/torus[512,k] for k in selected]
assert min(ratios)>3.9 and max(ratios)<4.2
assert torus_summary[-1]['max_error']<0.001
stored_report=None
if '--stored' in sys.argv:
path=Path(sys.argv[sys.argv.index('--stored')+1])
doc=json.loads(path.read_text())
records=doc['Numbers']
assert isinstance(records,list)
stored_keys=[(int(r['params']['m']),int(r['params']['n'])) for r in records]
assert stored_keys==selected
intervals={k:text_interval(str(r['number'])) for k,r in zip(stored_keys,records)}
with mp.workdps(150):
for k,(low,high) in intervals.items():
q=mp.mpf(independent[k,145])
if low==high:
assert k==(1,0) and low==QQ(1)/2
assert abs(q-mp_rational(low))<mp.mpf('1e-135')
else:
assert mp_rational(low)<q<mp_rational(high)
intervals[0,0]=(QQ(0),QQ(0))
count=0
for m,n in selected:
neighbors=[orbit(m+1,n),orbit(m-1,n),orbit(m,n+1),orbit(m,n-1)]
if not all(k in intervals for k in neighbors):continue
low,high=intervals[m,n]
residual_low=4*low-sum(intervals[k][1] for k in neighbors)
residual_high=4*high-sum(intervals[k][0] for k in neighbors)
assert residual_low<=0<=residual_high
count+=1
for k in range(1,11):
with mp.workdps(150):
expected=2/mp.pi*sum(mp.mpf(1)/(2*j-1) for j in range(1,k+1))
lo,hi=intervals[k,k]
assert mp_rational(lo)<expected<mp_rational(hi)
stored_report={'all_120_independent_integrals_inside_stored_intervals':True,
'interior_harmonic_identities':count,'diagonal_identities':10,
'exact_rational':'r(1,0)=1/2','stored_sha256':hashlib.sha256(path.read_bytes()).hexdigest()}
report={'table':'T596','manifest':'agents/runs/20261029T011133Z-build-skills/manifest.json',
'skills_sha256':'28ecbeeccb0c5e5d6eb541b5a5b0f88aec47b154374d909085116d34f5629199',
'versions':{'python':platform.python_version(),'sage':sage.version.version,
'numberdb':client_module.__version__,'mpmath':mp.__version__,'numpy':np.__version__},
'entries':120,'subsets':counts,'exact_entries':1,
'size':measure(written),'minimum_conservative_relative_bits':min_accuracy_bits,
'maximum_working_bits':max_precision,'integral_controls':controls,
'quadrature_agreement_absolute_tolerance':'1e-112',
'max_recurrence_integral_residual_145_digits':mp.nstr(max_residual,10),
'exact_harmonic_identities':exact_harmonic,'exact_diagonal_identities':10,
'torus_checks':torus_summary,'torus_256_to_512_error_ratio_range':[min(ratios),max(ratios)],
'stored_checks':stored_report,'rows':rows}
print(json.dumps(report,indent=2))
if __name__=='__main__':
run()