back to table · edit · history · where entries came from · files · download
9136 bytes, as of the version from 2026-10-25 11:31 (current). Recorded here, not run.
"""Independent checks for numberdb.org/T484.
sage -python checks.py
sage -python checks.py --stored # NUMBERDB_API_KEY required for draft
The printed estimates were not obtained by the finite enumerations here.
Enumeration checks the conventions; published independent estimates, bounds,
and a series extrapolation check the limiting values. None of these certifies
the Monte Carlo or series-analysis error estimates. Conversion is checked by
Python Decimal independently of the generator's Arb arithmetic.
"""
import os
import sys
import json
import itertools
from pathlib import Path
from fractions import Fraction
from decimal import Decimal, localcontext
import numberdb.sage as numberdb
from sage.rings.rational_field import QQ
from generate import AnimalGrowth, cubic_site_text
from numberdb._generate import _stored_entries
ROOT = Path(__file__).resolve().parent
DATA = json.loads((ROOT/'source-data.json').read_text())
def interval(text):
center, radius = map(Fraction, text.split(' +/- '))
assert radius > 0
return center-radius, center+radius
def report(label, value):
print(label+':', value, flush=True)
def canonical(points):
shift = tuple(min(p[i] for p in points) for i in range(len(next(iter(points)))))
return tuple(sorted(tuple(p[i]-shift[i] for i in range(len(shift))) for p in points))
def neighbours(point):
for i in range(len(point)):
for sign in [-1,1]:
q=list(point);q[i]+=sign
yield tuple(q)
def site_sets(d, maximum):
current={((0,)*d,)}
result={1:current}
for n in range(2,maximum+1):
following=set()
for animal in current:
points=set(animal)
boundary={y for x in animal for y in neighbours(x)}-points
for y in boundary: following.add(canonical(points|{y}))
current=following;result[n]=current
return result
def canonical_edges(edges):
points={p for edge in edges for p in edge}
d=len(next(iter(points)))
shift=tuple(min(p[i] for p in points) for i in range(d))
return tuple(sorted(tuple(sorted(tuple(p[i]-shift[i] for i in range(d)) for p in edge)) for edge in edges))
def bond_counts(d, maximum):
origin=(0,)*d
current={canonical_edges({(origin,q)}) for q in neighbours(origin)}
counts={1:len(current)}
for n in range(2,maximum+1):
following=set()
for animal in current:
edges=set(animal); vertices={p for e in animal for p in e}
for p in vertices:
for q in neighbours(p):
edge=tuple(sorted([p,q]))
if edge not in edges:
following.add(canonical_edges(edges|{edge}))
current=following; counts[n]=len(current)
return counts
def induced_subset_counts(site_animals, maximum):
"""Independent edge-counting algorithm: subsets of induced site graphs."""
counts={m:0 for m in range(1,maximum+1)}
for n,animals in site_animals.items():
if n<2 or n>maximum+1: continue
for animal in animals:
vertices=set(animal)
edges=sorted((p,q) for p in animal for q in neighbours(p) if p<q and q in vertices)
for m in range(n-1,min(maximum,len(edges))+1):
for selected in itertools.combinations(edges,m):
reached={animal[0]}
while True:
expanded=reached|{v for e in selected if any(u in reached for u in e) for v in e}
if expanded==reached: break
reached=expanded
if reached==vertices: counts[m]+=1
return counts
def extrapolate(counts, indices):
"""Extrapolate n*a_n/((n-1)*a_(n-1)) as polynomial in 1/n."""
result=QQ(0)
for n in indices:
value=QQ(n)*QQ(counts[n])/(QQ(n-1)*QQ(counts[n-1]))
x=QQ(1)/QQ(n)
weight=QQ(1)
for m in indices:
if m!=n:
y=QQ(1)/QQ(m);weight*=(-y)/(x-y)
result+=weight*value
return result
def records_from_document(document):
# Reuse the client's parser: the API normalizes flat submissions to a
# parameter-nested Numbers mapping. Check the exact identity set too.
numbers = _stored_entries(document)
g = AnimalGrowth()
expected = {','.join(p[k] for k in g.parameters) for p in g.enumerate()}
assert set(numbers) == expected, 'missing or extra stored entry'
return [dict(params=p, number=numbers[','.join(p[k] for k in g.parameters)])
for p in g.enumerate()]
def run(stored=None):
g=AnimalGrowth()
entries=stored if stored is not None else [dict(g.value(p,g.digits),params=p) for p in g.enumerate()]
assert [e['params'] for e in entries]==list(g.enumerate())
values={(e['params']['class'],e['params']['lattice']):interval(e['number']) for e in entries}
assert len(values)==4
# Independent arithmetic control first, including a known nonzero result.
with localcontext() as context:
context.prec=70
assert abs(Decimal(2).ln().exp()-Decimal(2))<Decimal('1e-65')
a=Decimal('2.1218588');r=Decimal('0.0000025')
low=(a-r).exp();high=(a+r).exp()
lo,hi=values['site','sc']
dlo=Decimal(lo.numerator)/Decimal(lo.denominator)
dhi=Decimal(hi.numerator)/Decimal(hi.denominator)
assert dlo<low<high<dhi
assert dhi-dlo < Decimal('1.1')*(high-low)
report('Decimal exponentiation control', 'exp(log(2)) agrees with 2 within 1e-65')
report('source log interval exponentiated independently',[str(low),str(high)])
a2=Decimal('1.4018155');r2=Decimal('0.0000030')
mc=(Fraction((a2-r2).exp()),Fraction((a2+r2).exp()))
sq=values['site','square']
assert mc[0]<sq[0]<sq[1]<mc[1]
old_sc=(Fraction(Decimal('2.118').exp()),Fraction(Decimal('2.122').exp()))
assert old_sc[0]<lo<hi<old_sc[1]
assert Fraction('2.1218563') < Fraction('2.1218572') < Fraction('2.1218612') < Fraction('2.1218613')
assert cubic_site_text(192)==cubic_site_text(384)
report('Arb precision repeat', '192 and 384 bits give identical written interval')
for key,ref in [(('bond','square'),('5.208','0.004')),(('bond','sc'),('10.62','0.08'))]:
lo,hi=values[key];c,r=map(Fraction,ref)
assert c-r<lo<hi<c+r
report('independent published estimates', 'all four selected intervals contained in earlier independent estimates; bond and cubic series values read as quoted by HNG')
for c,r in DATA['newer_square_estimates_not_selected']:
c,r=map(Fraction,(c,r));assert sq[0]<c-r<c+r<sq[1]
assert Fraction('4.06256912')-Fraction('0.00000002')>Fraction('4.06256903')+Fraction('0.00000005')
report('newer square estimates', 'non-overlap confirmed; neither selected, both contained in retained older estimate')
assert Fraction('4.00253')<sq[0]<sq[1]<=Fraction('4.5238')
assert values['site','sc'][1]<Fraction('9.3835')
for (cl,lat),(lo,hi) in values.items():
d=2 if lat=='square' else 3
assert d<lo<hi<(2*d)**2
report('published rigorous site bounds and elementary size bounds','passed')
# Controls for each exact enumerator have known answers.
assert [len(x) for x in site_sets(1,6).values()]==[1]*6
assert list(bond_counts(1,5).values())==[1]*5
report('enumeration controls','one-dimensional site and bond counts all equal 1')
for d,max_n,lat in [(2,8,'square'),(3,6,'sc')]:
animals=site_sets(d,max_n)
got={n:len(s) for n,s in animals.items()}
expected={n:int(DATA['counts'][lat+'-site'][str(n)]) for n in got}
assert got==expected
bc=bond_counts(d,4)
other=induced_subset_counts(animals,4)
assert bc==other
if d==2: assert list(bc.values())==[2,6,22,88]
report(lat+' exact site enumeration',got)
report(lat+' bond enumeration, both algorithms',bc)
counts={int(n):QQ(v) for n,v in DATA['counts']['square-site'].items()}
assert int(DATA['count_correction']['published'])-int(DATA['count_correction']['corrected'])==2
synthetic={n:QQ(7)**n/QQ(n) for n in range(1,71)}
indices=[40,45,50,55,60,65,70]
assert extrapolate(synthetic,indices)==7
report('series extrapolation control','a_n=7^n/n returns exactly 7')
for indices in [[30,35,40,45,50,55,60],[40,45,50,55,60,65,70]]:
estimate=extrapolate(counts,indices)
assert QQ(sq[0].numerator)/sq[0].denominator < estimate < QQ(sq[1].numerator)/sq[1].denominator
report('independent series extrapolation '+str(indices),str(estimate.numerical_approx(digits=18)))
report('precision limitation','extrapolation and source checks are heuristic; no final digits of any growth constant are certified')
report('RESULT','PASS')
if __name__=='__main__':
use_stored='--stored' in sys.argv or os.environ.get('CHECK_STORED')=='1'
if os.environ.get('NUMBERDB_KEY_FROM_STDIN')=='1':
numberdb.configure(api_key=sys.stdin.read().strip())
doc=numberdb.table('T484') if use_stored else None
run(records_from_document(doc) if doc else None)