back to table · edit · history · where entries came from · files · download
7431 bytes, as of the version from 2026-10-29 09:12 (current). Recorded here, not run.
"""Independent numerical checks for numberdb.org/T605.
Run: sage -python checks.py [--stored]
The optional stored mode reads the draft with NUMBERDB_API_KEY. It compares
all written claims, not just their midpoints, against a separately evaluated
four-sign average and a Gaussian integral. mpmath quadrature is a diagnostic,
not the source of the table's proven error bounds. Arb supplies those bounds.
Controls test Gaussian moments 0, 2, 4 and an exact four-sign moment before
using quadrature on selected exponents. The finite-vector inequality checks
are diagnostics, not a proof of Haagerup's theorem.
"""
import sys
import os
import json
import platform
from fractions import Fraction
from decimal import Decimal
from itertools import product
from math import comb
import numberdb.sage as numberdb
from numberdb._write import Entries
from sage.rings.rational_field import QQ
from sage.rings.real_arb import RealBallField
from sage.env import SAGE_VERSION
import mpmath as mp
from generate import KhintchineLower, evaluate, selection_groups, WORKING_GUARD
def to_mp(q):
q = Fraction(str(q))
return mp.mpf(q.numerator) / q.denominator
def claim(s):
if '.' not in s and 'e' not in s.lower():
q = Fraction(s)
return q, q, None
d = Decimal(s)
ulp = Fraction(10) ** d.as_tuple().exponent
q = Fraction(d)
return q-ulp, q+ulp, ulp
def gaussian_moment(p):
# x=t^2 regularizes the endpoint. No Gamma evaluation is used.
f = lambda t: 2*t**(2*p+1)*mp.exp(-t**4/2)
return mp.sqrt(2/mp.pi)*mp.quad(f, [0, 1, 2, 4, mp.inf])
def sign_norm(p, coefficients):
norm2 = mp.sqrt(sum(mp.mpf(a)**2 for a in coefficients))
moment = mp.fsum(abs(sum(mp.mpf(a)*e for a,e in zip(coefficients, signs)))**p
for signs in product((-1,1), repeat=len(coefficients)))
moment /= 2**len(coefficients)
return moment**(1/p)/norm2
def independent(p):
if p >= 2:
return mp.mpf(1), 'unit'
two = sign_norm(p, (1,1))
gaussian = gaussian_moment(p)**(1/p)
if two < gaussian:
return two, 'two-point'
return gaussian, 'gaussian'
def run(stored=False):
mp.mp.dps = 150
controls = {}
for p, known in [(0,1), (2,1), (4,3)]:
error = abs(gaussian_moment(mp.mpf(p))-known)
assert error < mp.mpf('1e-140'), (p,error)
controls['gaussian_moment_%s' % p] = mp.nstr(error,8)
assert abs(sign_norm(mp.mpf(1),(1,1))-1/mp.sqrt(2)) < mp.mpf('1e-140')
assert abs(sign_norm(mp.mpf(2),(1,2,3))-1) < mp.mpf('1e-140')
controls['sign_average'] = 'p=1 at (1,1) and p=2 at (1,2,3) passed'
# A ten-ulp mutation must fail the written-claim test.
original = mp.mpf(1)/mp.sqrt(2)
s = mp.nstr(original,100)
lo,hi,ulp=claim(s)
assert to_mp(lo) < original < to_mp(hi)
assert not (to_mp(lo+10*ulp) < original < to_mp(hi+10*ulp))
controls['ten_ulp_mutation_rejected'] = True
gen = KhintchineLower()
fields = Entries('p')
generated = {}
for params in gen.enumerate():
v = gen.value(params,100)
annotated = v if isinstance(v,dict) else {'number':v}
fields.add(p=params['p'], digits=100, **annotated)
generated[params['p']] = annotated['number']
records = fields.as_list()
if stored:
doc = numberdb.table('T605')
records = doc['Numbers']
if isinstance(records, dict):
records = [dict({'params': {'p': str(p)}}, **(v if isinstance(v, dict) else {'number': v}))
for p, v in records.items()]
assert isinstance(records, list)
assert [r['params']['p'] for r in records] == list(generated)
groups = {k:[str(x) for x in v] for k,v in selection_groups().items()}
assert len(set().union(*map(set,groups.values()))) == len(records) == 36
assert [Fraction(r['params']['p']) for r in records] == sorted(Fraction(p) for p in generated)
evidence=[]
previous=None
for rec in records:
ptext=rec['params']['p']; q=QQ(ptext); p=to_mp(ptext)
text=str(rec['number']); lo,hi,ulp=claim(text)
answer, independent_branch=independent(p)
value, branch=evaluate(q,100)
assert branch.startswith(independent_branch) or branch==independent_branch
assert mp.isfinite(answer) and answer>0
assert previous is None or answer>previous
previous=answer
if ulp is None:
assert value == QQ(text)
assert abs(answer-to_mp(text)) < mp.mpf('1e-140')
normalized_error='exact'
relative_radius='0'
radius='0'
else:
assert to_mp(lo)<answer<to_mp(hi), (ptext,text,mp.nstr(answer,110))
normalized_error=mp.nstr(abs(answer-(to_mp(lo)+to_mp(hi))/2)/to_mp(ulp),12)
field=value.parent()
assert value > field(QQ(lo.numerator)/QQ(lo.denominator))
assert value < field(QQ(hi.numerator)/QQ(hi.denominator))
relative=field(value.rad())/abs(value)
assert relative < field(10)**(-125), (ptext,relative)
relative_radius=str(relative)
radius=str(value.rad())
if branch.startswith('two-point'):
a,b=q.numerator(),q.denominator()
polynomial_power=RealBallField(512)(value)**(2*a)
target=QQ(2)**(a-2*b)
assert target in polynomial_power
# Separate, explicit vectors exercise the probability normalization.
for vector in [(1,), (1,1), (1,2), (1,1,1), (1,2,3)]:
assert sign_norm(p,vector)+mp.mpf('1e-140')>=answer
row={'p':ptext,'number':text,'branch':branch,
'radius':radius,'relative_radius':relative_radius,
'independent_error_in_last_places':normalized_error}
evidence.append(row)
print('checked p='+ptext, file=sys.stderr, flush=True)
for q in [QQ(2),QQ(5)/2,QQ(3),QQ(10)]:
assert evaluate(q)[0]==1
return {'table':'T605','stored_readback':stored,'entries':len(records),
'skill_manifest':'agents/runs/20261029T085909Z-build-skills/manifest.json',
'skill_sha256':'28ecbeeccb0c5e5d6eb541b5a5b0f88aec47b154374d909085116d34f5629199',
'working_bits':numberdb.bits(100,losing=WORKING_GUARD),
'relative_radius_bound':'< 1e-125 for every inexact computed value',
'quadrature_dps':150,'controls':controls,'groups':groups,
'versions':{'python':platform.python_version(),'sage':SAGE_VERSION,
'mpmath':mp.__version__,'numberdb':getattr(numberdb,'__version__','unknown')},
'records':records,'evidence':evidence,
'checks_passed':['exactness and named imports via dry_run.py',
'strict Arb branch comparisons', '100-digit written claims contain Arb balls',
'independent four-sign and Gaussian-integral evaluation for every row',
'algebraic power identity on every two-point row',
'five finite coefficient vectors at every p',
'strict monotonicity on selected inputs',
'exact unit-branch controls at p=2,5/2,3,10']}
if __name__=='__main__':
if os.environ.get('NUMBERDB_KEY_FROM_STDIN')=='1':
os.environ['NUMBERDB_API_KEY']=sys.stdin.read().strip()
print(json.dumps(run('--stored' in sys.argv),indent=2))