back to table · edit · history · where entries came from · files · download
12309 bytes, as of the version from 2026-10-27 12:26 (current). Recorded here, not run.
"""Independent numerical checks for Laplacian energies -- numberdb.org/T548.
Run: sage -python checks.py [--stored]
The default checks generated 100-digit spellings. --stored reads the draft
using NUMBERDB_API_KEY, or the public table after publication. No writes.
Place the generator, census-reference.json and corpus-reference.json beside
this file. All external inputs are archived snapshots with identified sources.
The graph6 decoder and permutation determinant are independent of Sage's
graph matrices. Exact rational Sturm isolation and bisection enclose every
root without AA, and sum absolute deviations directly. mpmath's symmetric
matrix eigensolver supplies a second algorithm at 150 decimal digits.
"""
import os
import sys
import json
import hashlib
from pathlib import Path
from itertools import permutations
from collections import Counter
from fractions import Fraction
from decimal import Decimal
from functools import lru_cache
import numberdb.sage as numberdb
from numberdb._write import to_text
from sage.rings.rational_field import QQ
from sage.rings.integer_ring import ZZ
from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
from sage.rings.real_mpfi import RealIntervalField
from sage.graphs.graph import Graph
from sage.env import SAGE_VERSION
import mpmath as mp
from generate import LaplacianEnergies, graph_keys, laplacian_energy, spectral_parts
HERE = Path(__file__).resolve().parent
R = PolynomialRing(QQ, 'x')
x = R.gen()
TARGET = QQ(1) / ZZ(2)**440
def qfrac(f):
return QQ(f.numerator) / QQ(f.denominator)
def written_bounds(text):
text = str(text)
c = qfrac(Fraction(text))
if '.' not in text and 'e' not in text.lower():
return c, c
e = Decimal(text).as_tuple().exponent
unit = QQ(10)**e
return c-unit, c+unit
def number_of(value):
return value['number'] if isinstance(value, dict) else value
def decode(key):
n = ord(key[0])-63
assert 0 <= n <= 62
bits = [((ord(c)-63) >> b) & 1 for c in key[1:] for b in range(5,-1,-1)]
A = [[0]*n for _ in range(n)]
k = 0
for j in range(1,n):
for i in range(j):
A[i][j] = A[j][i] = bits[k]
k += 1
assert not any(bits[k:])
return A
def laplacian(A):
return [[sum(A[i]) if i == j else -A[i][j] for j in range(len(A))]
for i in range(len(A))]
def determinant_polynomial(M):
n = len(M)
out = R(0)
for perm in permutations(range(n)):
sign = -1 if sum(perm[i]>perm[j] for i in range(n) for j in range(i+1,n))%2 else 1
term = R(sign)
for i,j in enumerate(perm):
term *= x-M[i][j] if i == j else -M[i][j]
if not term:
break
out += term
return out
def variations(sequence, point):
signs = [v.sign() for f in sequence if (v := f(point)) != 0]
return sum(a!=b for a,b in zip(signs, signs[1:]))
@lru_cache(maxsize=None)
def isolate(f, upper):
# f is irreducible. Linear factors are handled exactly.
if f.degree() == 1:
root = -f[0]/f[1]
return ((root,root),)
seq = [f, f.derivative()]
while True:
remainder = -(seq[-2] % seq[-1])
if not remainder:
break
seq.append(remainder)
left,right = QQ(0),QQ(upper)
vl,vr = variations(seq,left),variations(seq,right)
assert vl-vr == f.degree()
stack = [(left,right,vl,vr)]
intervals = []
while stack:
lo,hi,a,b = stack.pop()
if a-b == 0:
continue
if a-b == 1:
flo = f(lo)
assert flo*f(hi) < 0
while hi-lo > TARGET:
mid = (lo+hi)/QQ(2)
fm = f(mid)
assert fm != 0 # no rational roots of a nonlinear irreducible
if flo*fm < 0:
hi = mid
else:
lo,flo = mid,fm
assert variations(seq,lo)-variations(seq,hi) == 1
intervals.append((lo,hi))
else:
mid = (lo+hi)/QQ(2)
v = variations(seq,mid)
stack.extend(((lo,mid,a,v),(mid,hi,v,b)))
return tuple(sorted(intervals))
def independent_energy(p,n,m):
d = QQ(2*m)/QQ(n)
low=high=QQ(0)
root_moments = [QQ(0),QQ(0),QQ(0),QQ(0)]
all_rational = True
for f,mult in p.factor():
intervals = isolate(f,n+1)
assert len(intervals) == f.degree()
for a,b in intervals:
root_moments[0] += mult*a
root_moments[1] += mult*b
root_moments[2] += mult*a*a
root_moments[3] += mult*b*b
if all(a>=d for a,b in intervals):
exact = mult*(-f[f.degree()-1]/f.leading_coefficient()-f.degree()*d)
low+=exact; high+=exact
elif all(b<=d for a,b in intervals):
exact = mult*(f.degree()*d+f[f.degree()-1]/f.leading_coefficient())
low+=exact; high+=exact
else:
all_rational=False
for a,b in intervals:
assert b<d or a>d
u,v = (a-d,b-d) if a>d else (d-b,d-a)
low+=mult*u; high+=mult*v
assert root_moments[0] <= 2*m <= root_moments[1]
return low,high,all_rational,root_moments
def mpq(q):
return mp.mpf(str(q.numerator()))/mp.mpf(str(q.denominator()))
def diagonal_energy(A):
n=len(A);L=laplacian(A)
d=mp.mpf(sum(map(sum,A)))/n
eig=mp.eigsy(mp.matrix(L),eigvals_only=True)
assert all(mp.isfinite(v) for v in eig)
return sum(abs(v-d) for v in eig)
def main(stored=False):
mp.mp.dps=150
# Known-answer and deliberately false-answer controls run first.
assert decode('DQc') == [[0,0,1,0,1],[0,0,0,1,0],[1,0,0,0,0],[0,1,0,0,1],[1,0,0,1,0]]
assert determinant_polynomial(laplacian(decode('A_'))) == x*(x-2)
assert independent_energy(x*(x-2),2,1)[:3] == (QQ(2),QQ(2),True)
assert abs(diagonal_energy(decode('A_'))-2) < mp.mpf('1e-140')
assert diagonal_energy([[0,0],[0,0]]) == 0
repeated=(x-2)**3*x
assert sum(mult*f.degree() for f,mult in repeated.factor()) == 4
# A polynomial with two positive irrational roots tests Sturm isolation.
z=isolate(x*x-4*x+2,5)
assert len(z)==2 and z[0][1]<1 and 3<z[1][0]
I=RealIntervalField(413);third=I(QQ(1)/QQ(3))
assert third.lower().exact_rational()<third.upper().exact_rational()
assert written_bounds('2') == (QQ(2),QQ(2))
assert not (written_bounds('3')[0] <= 2 <= written_bounds('3')[1])
print('Known-answer controls passed.',flush=True)
keys=graph_keys()
corpus=json.loads((HERE/'corpus-reference.json').read_text())
census=json.loads((HERE/'census-reference.json').read_text())
all_catalogue=set()
for n,rawkeys in census['graphs'].items():
canonical={Graph(k,format='graph6').canonical_label(algorithm='sage').graph6_string() for k in rawkeys}
assert len(canonical)==len(rawkeys)
assert canonical=={k for k in keys if ord(k[0])-63==int(n)}
all_catalogue.update(canonical)
assert all_catalogue==set(keys)
for tid,doc in corpus.items():
assert set(doc['numbers'])==set(keys),tid
g=LaplacianEnergies()
raw={key:g.value({'g':key},100) for key in keys}
values={key:to_text(number_of(v),100) for key,v in raw.items()}
if stored:
doc=numberdb.table('T548')
values={key:str(number_of(v)) for key,v in doc['Numbers'].items()}
assert list(values)==list(keys)
report={'table':'T548','stored_readback':stored,'sage_version':SAGE_VERSION,
'mpmath_version':mp.__version__,'working_bits':413,'decimal_digits':100,
'skill_manifest':'agents/runs/20261027T121056Z-build-skills/manifest.json',
'skill_sha256':'e37efb7e2ededf6c1d7134c673f152f648f26e0267c556aebdf3ec65765d8b4c',
'counts':dict(Counter(ord(k[0])-63 for k in keys)),
'controls':'graph6 worked example; K2 spectrum and energy; edgeless matrix; Sturm isolator; rational endpoint width; wrong value rejected',
'exact':0,'irrational':0,'regular':0,'tree':0,'bipartite':0,
'formula_checks':{'complete':0,'path':0,'cycle':0,'complete_bipartite':0},
'rows':{}}
max_width=QQ(0);max_error=mp.mpf(0)
for key in keys:
A=decode(key);n=len(A);m=sum(map(sum,A))//2;deg=list(map(sum,A))
G=Graph(key,format='graph6')
assert G.canonical_label(algorithm='sage').graph6_string()==key
assert A==[list(row) for row in G.adjacency_matrix()]
L=laplacian(A);p=determinant_polynomial(L)
expected=R(str(number_of(corpus['T238']['numbers'][key])))
assert p==expected==R(G.laplacian_matrix().charpoly()),key
low,high,rational,moments=independent_energy(p,n,m)
assert moments[2] <= 2*m+sum(d*d for d in deg) <= moments[3]
a,b=written_bounds(values[key])
assert a<=low<=high<=b,(key,'independent enclosure outside stored spelling')
assert rational==(a==b),(key,'rationality classification')
if rational:
report['exact']+=1
else:
report['irrational']+=1
v=number_of(raw[key]);width=v.upper().exact_rational()-v.lower().exact_rational()
max_width=max(max_width,width)
assert width>0
assert a <= v.lower().exact_rational() <= v.upper().exact_rational() <= b
numerical=diagonal_energy(A)
err=abs(numerical-(mpq(low)+mpq(high))/2)
max_error=max(err,max_error)
assert err<mp.mpf('1e-125'),key
# The positive-tail relation follows from trace zero; validate it
# numerically independently of the root-isolation generator.
eig=mp.eigsy(mp.matrix(L),eigvals_only=True);d=mp.mpf(2*m)/n
assert abs(numerical-2*sum(v-d for v in eig if v>d))<mp.mpf('1e-140')
if len(set(deg))==1:
report['regular']+=1
ea,eb=written_bounds(number_of(corpus['T546']['numbers'][key]))
assert max(a,ea)<=min(b,eb)
adj=mp.eigsy(mp.matrix(A),eigvals_only=True)
assert abs(numerical-sum(abs(t) for t in adj))<mp.mpf('1e-140')
if m==n-1:report['tree']+=1
if G.is_bipartite():report['bipartite']+=1
if m==n*(n-1)//2:
assert a==b==2*(n-1);report['formula_checks']['complete']+=1
if m==n-1 and max(deg)<=2:
expected=sum(abs(2-2*mp.cos(mp.pi*j/n)-d) for j in range(n))
assert abs(numerical-expected)<mp.mpf('1e-140');report['formula_checks']['path']+=1
if n>=3 and all(v==2 for v in deg):
expected=sum(abs(2-2*mp.cos(2*mp.pi*j/n)-2) for j in range(n))
assert abs(numerical-expected)<mp.mpf('1e-140');report['formula_checks']['cycle']+=1
if G.is_bipartite():
s,t=sorted(map(len,G.bipartite_sets()))
if m==s*t:
mean=QQ(2*s*t)/QQ(n)
expected=n+(s-1)*abs(t-mean)+(t-1)*abs(s-mean)
assert a==b==expected
report['formula_checks']['complete_bipartite']+=1
report['rows'][key]={'number':values[key],'exact':rational,
'independent_lower':str(low),'independent_upper':str(high),
'mp_error':mp.nstr(err,8)}
report['max_generator_width']=str(max_width)
report['max_generator_width_decimal']=mp.nstr(mpq(max_width),12)
report['max_independent_diagonalization_error']=mp.nstr(max_error,12)
report['longest_written_value']=max(map(len,values.values()))
report['corpus_key_checks']=list(corpus)
report['all_142_permutation_determinants_match_T238']=True
report['all_142_sturm_enclosures_within_written_intervals']=True
report['all_85_generator_enclosures_within_written_intervals']=True
report['all_142_mpmath_diagonalizations_agree']=True
report['all_142_trace_positive_tail_and_second_moments_pass']=True
print('All checks passed.',flush=True)
print('EVIDENCE_JSON '+json.dumps(report,sort_keys=True),flush=True)
return report
if __name__=='__main__':
if os.environ.get('NUMBERDB_KEY_FROM_STDIN')=='1':
key=sys.stdin.read().strip()
if '=' in key and key.split('=',1)[0].strip().isidentifier():key=key.split('=',1)[1].strip().strip('\"\'')
os.environ['NUMBERDB_API_KEY']=key
main(stored='--stored' in sys.argv)