back to table · edit · history · where entries came from · files · download
10834 bytes, as of the version from 2026-10-27 11:38 (current). Recorded here, not run.
"""Independent checks for Graph energies of connected graphs, numberdb.org/T546.
Run with all files alongside this script:
sage -python checks.py
sage -python checks.py --stored
The independent numerical path decodes graph6 without Sage's Graph decoder,
then uses mpmath's symmetric-matrix eigensolver at 150 and 180 decimal digits.
It never calls the generator's characteristic-polynomial or root routines.
It is an agreement check, not a second certification of the enclosures.
Exact integer walk traces independently verify characteristic polynomials.
McKay's five catalogue files and NumberDB's related tables are read online if
not saved beside the script. Draft reads require NUMBERDB_API_KEY.
The JSON output contains entries, measurements, per-row evidence and versions.
"""
import sys
import os
import json
import pathlib
import urllib.request
from fractions import Fraction
from collections import Counter
import numberdb.sage as numberdb
from numberdb._write import to_text
from sage.rings.integer_ring import ZZ
from sage.rings.rational_field import QQ
from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
from sage.rings.real_mpfi import RealIntervalField
from sage.graphs.graph import Graph
from sage.graphs.graph_generators import graphs
from sage.env import SAGE_VERSION
import mpmath as mp
from generate import GraphEnergies, energy, spectral_factors
BASE = pathlib.Path(__file__).resolve().parent
POLYS = PolynomialRing(ZZ, 'x')
COUNTS = {2:1, 3:2, 4:6, 5:21, 6:112}
def adjacency(key):
"""Read the short-order graph6 specification directly."""
n = ord(key[0])-63
assert 0 <= n <= 62
bits = ''.join(format(ord(c)-63,'06b') for c in key[1:])
assert len(bits) >= n*(n-1)//2
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] = int(bits[k]); k += 1
return a
def multiply(a,b):
n=len(a)
return [[sum(a[i][k]*b[k][j] for k in range(n)) for j in range(n)] for i in range(n)]
def trace_polynomial(a):
"""Newton identities from exact integer closed-walk counts."""
n=len(a); power=[[int(i==j) for j in range(n)] for i in range(n)]
traces=[0]
for k in range(1,n+1):
power=multiply(power,a); traces.append(sum(power[i][i] for i in range(n)))
coeff=[1]
for k in range(1,n+1):
numerator=-sum(coeff[k-i]*traces[i] for i in range(1,k+1))
assert numerator % k == 0
coeff.append(numerator//k)
return POLYS(list(reversed(coeff)))
def independent(a,dps):
with mp.workdps(dps):
spectrum=mp.eigsy(mp.matrix(a),eigvals_only=True)
value=sum(abs(spectrum[j]) for j in range(len(a)))
assert mp.isfinite(value)
return value, list(spectrum)
def fraction(q):
return Fraction(int(q.numerator()),int(q.denominator()))
def enclosure(value):
if value in ZZ:
# Membership of an interval in ZZ is not an exactness test.
# Only the actual Sage Integer class may take this branch.
if type(value).__name__=='Integer':
q=Fraction(int(value)); return q,q
lo,hi=value.endpoints()
return fraction(lo.exact_rational()), fraction(hi.exact_rational())
def written_interval(text):
text=str(text).rstrip('?')
centre=Fraction(text)
if '.' not in text and 'e' not in text.lower():
return centre,centre
mantissa,sep,exponent=text.lower().partition('e')
places=len(mantissa.partition('.')[2])
unit=Fraction(10)**((int(exponent) if sep else 0)-places)
return centre-unit,centre+unit
def mpf(q):
return mp.mpf(q.numerator)/q.denominator
def document(tid):
path=BASE/(tid+'.json')
return json.loads(path.read_text()) if path.exists() else numberdb.table(tid)
def numbers(doc):
block=doc['Numbers']
if isinstance(block,dict):
return block
return {row['params']['g']:{k:v for k,v in row.items() if k!='params'} for row in block}
def number(row):
return row['number'] if isinstance(row,dict) else row
def run(stored=False):
# Controls precede any measurements.
known=adjacency('DQc')
assert [(i,j) for j in range(1,5) for i in range(j) if known[i][j]]==[(0,2),(1,3),(0,4),(3,4)]
for a, expected in [([[0,1],[1,0]],2),([[0,0],[0,0]],0),([[0,1,1],[1,0,1],[1,1,0]],4)]:
v,_=independent(a,150)
with mp.workdps(190): assert abs(v-expected)<mp.mpf('1e-140')
assert trace_polynomial([[0,1],[1,0]])==POLYS([-1,0,1])
gen=GraphEnergies(); params=list(gen.enumerate()); keys=[p['g'] for p in params]
assert len(keys)==142 and len(keys)==len(set(keys))
assert keys==sorted(keys,key=lambda k:(ord(k[0])-63,k))
assert Counter(ord(k[0])-63 for k in keys)==COUNTS
source_keys=[]
for n,count in COUNTS.items():
filename='mckay-%d.g6'%n
p=BASE/filename
raw=p.read_bytes() if p.exists() else urllib.request.urlopen('https://users.cecs.anu.edu.au/~bdm/data/graph%dc.g6'%n,timeout=30).read()
lines=raw.decode('ascii').splitlines(); assert len(lines)==count
source_keys.extend(Graph(k,format='graph6').canonical_label(algorithm='sage').graph6_string() for k in lines)
assert set(source_keys)==set(keys) and len(source_keys)==len(keys)
siblings={tid:numbers(document(tid)) for tid in ['T237','T238','T397','T400','T545']}
for tid,rows in siblings.items():
subset={k for k in rows if 2<=ord(k[0])-63<=6}
assert subset==set(keys),(tid,'census mismatch')
for k in keys:
assert Graph(k,format='graph6').canonical_label(algorithm='sage').graph6_string()==k
stored_rows=numbers(numberdb.table('T546')) if stored else None
if stored: assert list(stored_rows)==keys
records=[]; evidence=[]; exact=[]; widths=[]; work_bits=[]; counts=Counter(); max_agreement=mp.mpf(0)
for p in params:
key=p['g']; a=adjacency(key); graph=Graph(key,format='graph6'); n=len(a); edges=sum(map(sum,a))//2
assert a==[[int(graph.has_edge(i,j)) for j in range(n)] for i in range(n)]
assert graph.is_connected() and not graph.has_loops() and not graph.has_multiple_edges()
cp=trace_polynomial(a)
assert cp==POLYS(number(siblings['T237'][key]).replace('^','**'))
assert cp==graph.adjacency_matrix().charpoly()
value=gen.value(p,100); v=value['number']; lo,hi=enclosure(v)
text=to_text(v,100)
if stored: text=str(number(stored_rows[key]))
wlo,whi=written_interval(text)
assert wlo<=lo<=hi<=whi,(key,'certified enclosure outside written interval')
if lo==hi:
assert '.' not in text and lo.denominator==1 and lo.numerator%2==0
exact.append(key)
else:
widths.append(hi-lo);work_bits.append(int(v.parent().precision()))
assert 0<hi-lo<Fraction(1,10**110)
assert sum(c.isdigit() for c in text.split('e')[0])==100
low,eigs150=independent(a,150); high,eigs=independent(a,180)
with mp.workdps(190):
difference=abs(low-high); max_agreement=max(max_agreement,difference)
assert difference<mp.mpf('1e-140')
assert abs(high-(mpf(lo)+mpf(hi))/2)<mp.mpf('1e-120')
if lo!=hi: assert mpf(wlo)<high<mpf(whi)
else: assert abs(high-mpf(lo))<mp.mpf('1e-170')
assert abs(sum(eigs))<mp.mpf('1e-170')
assert abs(sum(x*x for x in eigs)-2*edges)<mp.mpf('1e-170')
assert abs(high-2*sum(x for x in eigs if x>0))<mp.mpf('1e-170')
assert high+mp.mpf('1e-170')>=2*mp.sqrt(edges)
assert high-mp.mpf('1e-170')<=mp.sqrt(2*n*edges)
if edges==n*(n-1)//2:
assert lo==hi==2*(n-1);counts['complete']+=1
if graph.is_bipartite():
counts['bipartite']+=1
left,right=graph.bipartite_sets()
if edges==len(left)*len(right):
expected=2*mp.sqrt(len(left)*len(right))
assert abs(high-expected)<mp.mpf('1e-170');counts['complete_bipartite']+=1
degrees=list(map(sum,a))
if edges==n-1:
counts['trees']+=1
if max(degrees)<=2:
expected=sum(abs(2*mp.cos(mp.pi*j/(n+1))) for j in range(1,n+1))
assert abs(high-expected)<mp.mpf('1e-170');counts['paths']+=1
if all(d==2 for d in degrees):
expected=sum(abs(2*mp.cos(2*mp.pi*j/n)) for j in range(n))
assert abs(high-expected)<mp.mpf('1e-170');counts['cycles']+=1
if len(set(degrees))==1: counts['regular']+=1
ev={'g':key,'characteristic_polynomial':str(cp),'written':text,'exact':lo==hi,
'certified_lower':str(lo),'certified_upper':str(hi),
'independent_180_digits':mp.nstr(high,145),'independent_difference':mp.nstr(difference,8)}
evidence.append(ev)
record={'params':p,'number':text}
if value.get('comment'):record['comment']=value['comment']
if stored and value.get('comment'):assert stored_rows[key].get('comment')==value['comment']
records.append(record)
if stored:
report=gen.verify(sample=None)
assert report.ok,str(report)
verify=str(report)
else:verify=None
summary={'table':'T546','mode':'stored' if stored else 'pre-fill','entries':len(records),
'counts_by_order':COUNTS,'exact_count':len(exact),'exact_keys':exact,
'irrational_count':len(widths),'maximum_interval_width':str(max(widths)),
'working_bits':sorted(set(work_bits)), 'max_150_vs_180_difference':mp.nstr(max_agreement,8),
'controls':'graph6 DQc; zero matrix; K2 energy 2; K3 energy 4; K2 trace polynomial x^2-1',
'structural_counts':dict(counts),'sibling_key_checks':list(siblings),
'catalogue_set_equality':True,'polynomial_checks':len(records),
'all_stored_enclosures_contain_certified_intervals':True if stored else None,
'longest_written_value':max(len(r['number']) for r in records),
'entries_block_bytes':len(json.dumps(records).encode()),
'verify':verify,'versions':{'sage':SAGE_VERSION,'numberdb':numberdb.__version__,'mpmath':mp.__version__},
'skill_manifest':'agents/runs/20261027T111010Z-build-skills/manifest.json',
'skill_local_sha256':'e37efb7e2ededf6c1d7134c673f152f648f26e0267c556aebdf3ec65765d8b4c',
'skill_live_sha256':'e37efb7e2ededf6c1d7134c673f152f648f26e0267c556aebdf3ec65765d8b4c'}
print(json.dumps({'summary':summary,'records':records,'evidence':evidence},indent=2))
if __name__=='__main__':
if os.environ.get('NUMBERDB_KEY_FROM_STDIN') == '1':
numberdb.configure(api_key=sys.stdin.read().strip())
run('--stored' in sys.argv or bool(os.environ.get('CHECK_STORED')))