back to table · edit · history · where entries came from · files · download
9869 bytes, as of the version from 2026-10-25 09:13 (current). Recorded here, not run.
"""Independent numerical checks for numberdb.org/T480.
sage -python checks.py
NUMBERDB_STORED=1 sage -python checks.py # authenticated draft read
Uses scipy's convex hull only to select facets; plane distances are
recomputed by explicit vector products in real ball arithmetic at 192 and
320 bits. The source coordinates are parsed exactly and radially normalised.
The finite-precision hull topology is not represented as a proof of global
optimality. Synthetic polyhedral controls use known answers and a deliberately
perturbed tetrahedron. Geometry uses no radius function from generate.py.
"""
import numberdb.sage as numberdb
from sage.rings.real_arb import RealBallField
from sage.rings.rational_field import QQ
from sage.rings.real_mpfi import RealIntervalField
from numberdb._write import to_text
import collections
import importlib.metadata
import itertools
import json
import math
import os
import pathlib
import sys
import numpy as np
import scipy
from scipy.spatial import ConvexHull
from generate import SphericalCoverings, SOURCE, EXACT, cosine, source_angle
from check import exactness, measure, names_its_rings, agrees_with
P = pathlib.Path(__file__).resolve().parent
def dec(text):
from decimal import Decimal
x = Decimal(text)
sign, digits, exponent = x.as_tuple()
q = QQ(int(''.join(map(str, digits)))) * QQ(10)**exponent
return -q if sign else q
def written(text, B):
from decimal import Decimal
q = dec(text)
return B(q, QQ(10)**Decimal(text).as_tuple().exponent) if '.' in text or 'e' in text.lower() else B(q)
def encloses(outer, inner):
return outer.lower() <= inner.lower() and inner.upper() <= outer.upper()
def dot(a, b):
return sum(x*y for x,y in zip(a,b))
def cross(a, b):
return [a[1]*b[2]-a[2]*b[1], a[2]*b[0]-a[0]*b[2], a[0]*b[1]-a[1]*b[0]]
def normalize(a):
length = dot(a, a).sqrt()
return [x/length for x in a]
def geometry(raw, bits):
B = RealBallField(bits)
points = [normalize([B(dec(x)) if isinstance(x,str) else B(x) for x in row]) for row in raw]
floats = np.array([[float(x) for x in row] for row in points])
hull = ConvexHull(floats)
assert np.all(hull.equations[:,3] < 0)
assert len(hull.simplices) == 2*len(raw)-4
edge_counts = collections.Counter(tuple(sorted((int(a),int(b)))) for face in hull.simplices for a,b in itertools.combinations(face,2))
assert all(v == 2 for v in edge_counts.values())
assert len(raw)-len(edge_counts)+len(hull.simplices) == 2
distances = []
worst_violation = 0.0
for ids in hull.simplices:
a,b,c = [points[i] for i in ids]
u = cross([b[k]-a[k] for k in range(3)], [c[k]-a[k] for k in range(3)])
u = normalize(u)
h = dot(u,a)
if h < 0:
u = [-x for x in u]; h = -h
assert h > 0
# Check every vertex lies behind each proposed supporting plane.
violation = max(float((dot(u,v)-h).upper()) for v in points)
assert violation < 1e-12
worst_violation = max(worst_violation,violation)
distances.append(h)
lo = min(d.lower() for d in distances); hi = min(d.upper() for d in distances)
distance = B(RealIntervalField(bits)(lo, hi))
rho = distance.arccos()*180/B.pi()
chord = (2-2*distance).sqrt()
assert rho.is_finite() and chord.is_finite()
assert float(rho.rad()) < 1e-40
float_radius = np.degrees(np.arccos(min(-hull.equations[:,3])))
assert abs(float(rho)-float_radius) < 1e-10
return rho, chord, worst_violation
def synthetic(n, B):
if n == 4:
return [[1,1,1],[1,-1,-1],[-1,1,-1],[-1,-1,1]]
if n == 6:
return [[s if j==k else 0 for j in range(3)] for k in range(3) for s in (-1,1)]
if n in (5,7):
m=n-2
return [[0,0,1],[0,0,-1]]+[[ (2*B.pi()*k/m).cos(), (2*B.pi()*k/m).sin(), B(0)] for k in range(m)]
if n == 12:
phi=(1+B(5).sqrt())/2
return [[0,s,t*phi] for s in (-1,1) for t in (-1,1)]+[[t*phi,0,s] for s in (-1,1) for t in (-1,1)]+[[s,t*phi,0] for s in (-1,1) for t in (-1,1)]
m=(n-2)//2
z=1/(1+1/(B.pi()/m).cos());a=(1-z*z).sqrt()
return [[0,0,1],[0,0,-1]]+[[a*(B.pi()*(2*k+phase)/m).cos(),a*(B.pi()*(2*k+phase)/m).sin(),s*z] for s,phase in [(1,0),(-1,1)] for k in range(m)]
def controls():
B=RealBallField(640)
for n in sorted(EXACT):
a,c,_=geometry(synthetic(n,B),512)
# Independently expressed closed forms, including each entry's formula.
cs={4:B(1)/3,5:1/B(5).sqrt(),6:1/B(3).sqrt(),
7:((3+B(5).sqrt())/(11+B(5).sqrt())).sqrt(),
10:((1+B(2).sqrt())/(3+B(2).sqrt())).sqrt(),
12:((5+2*B(5).sqrt())/15).sqrt(),
14:((3+3*B(3).sqrt())/(7+3*B(3).sqrt())).sqrt()}
expected=cs[n].arccos()*180/B.pi()
assert a.overlaps(expected),(n,a,expected)
assert abs(float(a)-float(SOURCE['rows'][str(n)]['angle'])) < 1e-7
assert cosine(n,B).overlaps(cs[n])
tetra=synthetic(4,B);tetra[0][0]=QQ(6)/5
bad,_,_=geometry(tetra,192)
assert abs(float(bad)-float((B(1)/3).arccos()*180/B.pi())) > 0.1
# The reusable equality helper checks the known octahedral centre norms.
assert not agrees_with(dict(enumerate([sum(x*x for x in a) for a in synthetic(6,B)])),lambda key:1)
print('Known-answer controls: seven exact configurations and a perturbed tetrahedron passed.',flush=True)
def run(stored=False):
controls()
generator=SphericalCoverings();B=RealBallField(512)
values={};entries={}
if stored:
if os.environ.get('NUMBERDB_KEY_FROM_STDIN') == '1':
os.environ['NUMBERDB_API_KEY']=sys.stdin.read().strip()
tree=numberdb.table('T480')
numbers=tree['Numbers']
if isinstance(numbers,list):
entries={(int(e['params']['n']),e['params']['normalisation']):e['number'] for e in numbers}
else:
entries={(int(n),form):(entry['number'] if isinstance(entry,dict) else entry) for n,forms in numbers.items() for form,entry in forms.items()}
def check_text(value):
if isinstance(value,str):
assert chr(92)*2 not in value, 'Unexpected double LaTeX escapes'
elif isinstance(value,dict):
for v in value.values():check_text(v)
elif isinstance(value,list):
for v in value:check_text(v)
for section in ('Definition','Parameters','Comments','Formulas','Numbers'):
check_text(tree.get(section))
code=tree['Programs']['program-sage']['code']
compile(code,'table-program','exec')
assert len(entries)==260
else:
for params in generator.enumerate():
key=(int(params['n']),params['normalisation']);v=generator.value(params,generator.digits_for(params));values[key]=v
entries[key]=to_text(v['number'],v['digits'])
assert not exactness(values)
assert not names_its_rings(P/'generate.py')
evidence={'mode':'database' if stored else 'local','entries':len(entries),'configuration_count':127,'versions':{'numberdb':numberdb.__version__,'scipy':scipy.__version__,'numpy':np.__version__},'rows':{}}
min_accuracy=10000; max_exact_radius=0.0;prior_angle=None;prior_chord=None
for n in range(4,131):
raw=SOURCE['rows'][str(n)]['coordinates'];assert len(raw)==n
a192,c192,v192=geometry(raw,192);a320,c320,v320=geometry(raw,320)
assert a192.overlaps(a320) and c192.overlaps(c320)
printed=B(dec(SOURCE['rows'][str(n)]['angle']));ulp=B(QQ(1)/10**7)
assert printed-ulp < a320 < printed,(n,printed,a320)
angle=written(entries[n,'angle'],B);chord=written(entries[n,'chord'],B)
for v in (angle,chord):assert v.is_finite()
assert angle.overlaps(360*(chord/2).arcsin()/B.pi())
assert chord.overlaps(2*(angle*B.pi()/360).sin())
# Area of n caps must be at least the sphere's area.
assert n*(1-(angle*B.pi()/180).cos()) > 2
if prior_angle is not None:
assert angle < prior_angle and chord < prior_chord
prior_angle,prior_chord=angle,chord
if n in EXACT:
params={'n':n,'normalisation':'angle'}
for norm,stored_value in [('angle',angle),('chord',chord)]:
params['normalisation']=norm
low=generator.value(params,100)['number'];high=generator.value(params,160)['number']
assert low.overlaps(high) and encloses(stored_value,high)
min_accuracy=min(min_accuracy,int(math.floor(-math.log2(float(low.rad())/abs(float(low))))))
max_exact_radius=max(max_exact_radius,float(low.rad()))
assert abs(float(a320-angle)) < 1e-7
assert abs(float(c320-chord)) < 2e-9
else:
assert encloses(angle,a320), (n,'angle misses coordinates')
assert encloses(chord,c320), (n,'chord misses coordinates')
converted=2*(source_angle(n,B)*B.pi()/360).sin()
assert encloses(chord,converted),(n,'chord misses source range')
evidence['rows'][str(n)]={'coordinate_angle':to_text(a320,60),'coordinate_chord':to_text(c320,60),'printed_minus_coordinates':float(printed-a320),'facet_violation':max(v192,v320),'angle':entries[n,'angle'],'chord':entries[n,'chord']}
if n%25==0:print('Checked configurations through n=%d'%n,flush=True)
evidence['minimum_exact_accuracy_bits']=min_accuracy
evidence['maximum_exact_ball_radius']=max_exact_radius
evidence['max_source_discrepancy_degrees']=max(r['printed_minus_coordinates'] for r in evidence['rows'].values())
if values:evidence['measurement']=measure(values)
print('EVIDENCE_JSON',json.dumps(evidence,separators=(',',':')))
return evidence
if __name__=='__main__':
run(stored=os.environ.get('NUMBERDB_STORED')=='1' or '--stored' in sys.argv)