back to table · edit · history · where entries came from · files · download
7568 bytes, as of the version from 2026-10-29 19:59 (current). Recorded here, not run.
"""Independent numerical and stored-interval checks for numberdb.org/T611.
sage -python checks.py
sage -python checks.py --stored
Requires numberdb and mpmath. --stored requires a key while T611 is private.
The generator uses Arb Li_2; this file uses mpmath quadrature, without Li_2.
Antiprisms: integrate the defining log-sine integral independently.
Lobell: use the geometrically derived arccosh integral in Mednykh--Vesnin,
Lobell manifolds revised (2007), Theorem 2, divided by eight for a polyhedron.
This tests the volume formula, not only its special-function evaluation.
A sign discrepancy is recorded deliberately: equation (2) in that paper has
+Lambda(2*theta-pi/2), whereas Inoue Theorem 4.1 has minus. The former sign
contradicts the same paper's independently derived Theorem 2 integral and
known dodecahedral volume. We use Inoue's formula and verify the discrepancy.
Numerical quadrature is a cross-check, not the proof of 100 certified digits.
Controls run first: integral 1, log-sine at pi/4 vs Catalan, and independent
arccosh integral at alpha=pi/4 vs Catalan/2. All comparisons require finite
answers and stated tolerances. The same checks apply to every database row.
"""
import json
import os
import sys
from decimal import Decimal
import numberdb.sage as numberdb
from numberdb._write import to_text
from sage.rings.real_arb import RealBallField
from sage.rings.rational_field import QQ
import mpmath as mp
from generate import RightAngledVolumes, volume
def log_sine(x):
# Subtract the endpoint logarithm. The ratio has the continuous value 1.
def f(u):
if not u:
return mp.mpf(0)
return mp.log(mp.sin(x*u)/(x*u))
return x*(1-mp.log(2*x)) - x*mp.quad(f,[0,1])
def schlafli_integral(alpha):
# Split at the logarithmic singularity t=pi/4.
# Its explicit primitive prevents lost digits near the singular endpoint.
c=mp.pi/4
b=mp.pi/3
endpoint_tolerance=100*mp.eps
def h(t):
u=t-c
ratio=1 if not u else mp.sin(2*u)/(2*u)
radicand=mp.cos(t)**2-mp.cos(2*t)**2
# The integrand is real on [alpha,b]. mpmath can have an endpoint
# imaginary roundoff of order sqrt(eps); take the real part there.
assert radicand >= -endpoint_tolerance
return mp.log(mp.cos(t)+mp.sqrt(max(mp.mpf(0),radicand)))-mp.log(ratio)
def primitive(u):
return mp.mpf(0) if not u else u*(mp.log(2*abs(u))-1)
return mp.quad(h,[alpha,c,b]) - primitive(b-c)+primitive(alpha-c)
def independent(family,n):
if family=='antiprism':
return 2*n*(log_sine(mp.pi/4+mp.pi/(2*n))
+log_sine(mp.pi/4-mp.pi/(2*n)))
return n*schlafli_integral(mp.pi/n)
def decimal_interval(text):
value=Decimal(text)
ulp=Decimal(1).scaleb(value.as_tuple().exponent)
return mp.mpf(str(value)),mp.mpf(str(ulp))
def rational_decimal(text):
value=Decimal(text)
sign,digits,exponent=value.as_tuple()
integer=int(''.join(map(str,digits)))*(-1 if sign else 1)
return QQ(integer)*QQ(10)**exponent
def check_decimal(text,ball,independent_value):
centre,ulp=decimal_interval(text)
assert mp.isfinite(independent_value)
error=abs(centre-independent_value)
assert error < ulp, (text,mp.nstr(error/ulp,10))
# Check exact decimal endpoints contain the certified computation.
exact_centre=rational_decimal(text)
exact_ulp=QQ(10)**Decimal(text).as_tuple().exponent
real=RealBallField(600)
assert real(exact_centre-exact_ulp) < ball < real(exact_centre+exact_ulp)
return error/ulp
def stored_rows():
document=numberdb.table('T611')
block=document['Numbers']
if isinstance(block,list):
return [(str(r['params']['family']),int(r['params']['n']),r['number']) for r in block]
rows=[]
for family,entries in block.items():
entries=entries.get('numbers',entries)
for n,entry in entries.items():
if n=='param-latex':continue
rows.append((family,int(n),entry['number'] if isinstance(entry,dict) else entry))
return rows
def run(use_stored=False):
mp.mp.dps=135
tol=mp.mpf('1e-120')
controls={
'quadrature_constant':abs(mp.quad(lambda x:mp.mpf(1),[0,1])-1),
'log_sine_Catalan':abs(log_sine(mp.pi/4)-mp.catalan/2),
'Schlafli_ideal_control':abs(schlafli_integral(mp.pi/4)-mp.catalan/2)}
assert all(mp.isfinite(x) and x < tol for x in controls.values()),controls
generator=RightAngledVolumes()
parameters=[(p['family'],p['n']) for p in generator.enumerate()]
assert parameters==([('antiprism',n) for n in range(3,31)]
+[('lobell',n) for n in range(5,31)])
balls={(family,n):volume(family,n) for family,n in parameters}
texts={(family,n):to_text(b, digits=100) for (family,n),b in balls.items()}
if use_stored:
rows=stored_rows()
assert [(f,n) for f,n,t in rows]==parameters
texts={(f,n):t for f,n,t in rows}
output=[]
max_ulp=mp.mpf(0)
min_digits=1000.0
for family,n in parameters:
ball=balls[family,n]
expected=independent(family,n)
text=texts[family,n]
ulp_error=check_decimal(text,ball,expected)
max_ulp=max(max_ulp,ulp_error)
relative_radius=float(ball.rad())/abs(float(ball.mid()))
import math
min_digits=min(min_digits,-math.log10(relative_radius))
# Independent quadrature at a second precision detects instability.
with mp.workdps(115):
second=independent(family,n)
assert abs(expected-second) < mp.mpf('1e-105')*abs(expected)
output.append({'family':family,'n':n,'number':text,
'quadrature':mp.nstr(expected,120),
'error_in_last_place':mp.nstr(ulp_error,12)})
assert abs(independent('antiprism',3)-4*mp.catalan) < tol
assert abs(independent('antiprism',4)-independent('lobell',6)) < tol
# Independently published 100 decimal digits of L_5 (OEIS A321209).
known=mp.mpf('4.306207600730808652919837159842209449241185235412934114224949117819904343308238335156704866571730974')
assert abs(independent('lobell',5)-known) < mp.mpf('1e-92')
# Demonstrate the printed sign error in Mednykh--Vesnin equation (2).
n=5
theta=mp.asin(1/(2*mp.cos(mp.pi/n)))
incorrect=mp.mpf(n)/2*(2*log_sine(theta)+log_sine(theta+mp.pi/n)
+log_sine(theta-mp.pi/n)-log_sine(mp.pi/2-2*theta))
correct=independent('lobell',5)
assert abs(correct-incorrect)>1
return {'table':'T611','stored':use_stored,'entries':len(output),
'controls':{k:mp.nstr(v,8) for k,v in controls.items()},
'min_relative_radius_decimal_digits':min_digits,
'max_error_in_last_place':mp.nstr(max_ulp,12),
'quadrature_precisions':[135,115],
'second_precision_relative_tolerance':'1e-105',
'source_sign_check':{'Inoue_and_integral':mp.nstr(correct,30),
'Mednykh_Vesnin_equation_2_printed_sign':mp.nstr(incorrect,30)},
'values':output,
'skill_manifest':'agents/runs/20261029T183732Z-build-skills/manifest.json',
'skill_sha256_local_and_live':'28ecbeeccb0c5e5d6eb541b5a5b0f88aec47b154374d909085116d34f5629199',
'mpmath_version':mp.__version__}
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))