back to table · edit · history · where entries came from · files · download
8690 bytes, as of the version from 2026-10-24 22:26 (current). Recorded here, not run.
"""Independent checks of T470. Run: sage -python checks.py.
With NUMBERDB_API_KEY (or NUMBERDB_KEY_FROM_STDIN=1), also check all stored
entries and the n=3 relation with T468. The exact polygon integration below
uses arbitrary convex polygon vertices and shares no formula or helper with
the generator. Disk quadrature integrates signed hull-edge areas, rather than
counting vertices and subtracting missing area. It is diagnostic; Arb in the
generator supplies the certified bounds. Simulations use a fixed random seed.
"""
import json
import os
import sys
import numberdb.sage as numberdb
from sage.rings.rational_field import QQ
from sage.rings.real_arb import RealBallField
from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
from check import exactness, measure, agrees_with, names_its_rings
from generate import ExpectedHullVolume, expected_ratio, REGIONS, MAX_N
import mpmath as mp
def cross(a, b):
return a[0]*b[1] - a[1]*b[0]
def polygon_ratio(vertices, n):
"""Integrate the signed area contributed by every possible hull edge.
Boundary chord A(s)->B(t) intersects two different polygon sides. If
D=B-A, the endpoint change of variables has factor
-cross(D,A')*cross(D,B')*(v-u) for 0<u<v<1. The signed cross product of
the two sampled points is (v-u)*cross(A,B); their remaining n-2 points
are on the left with probability (left_area/area)^(n-2).
Integral of (v-u)^2 on the unit triangle is 1/12. The shoelace formula
then supplies n(n-1)/2. Every other integration is an exact polynomial
integral over 0<s,t<1.
"""
P = PolynomialRing(QQ, names=('s', 't'))
s,t = P.gens()
vertices = [tuple(map(QQ,v)) for v in vertices]
m = len(vertices)
area = sum(cross(vertices[i],vertices[(i+1)%m]) for i in range(m))/2
assert area > 0
total = QQ(0)
for i in range(m):
for j in range(m):
if i == j:
continue
a0, a1 = vertices[i], vertices[(i+1)%m]
b0, b1 = vertices[j], vertices[(j+1)%m]
da = (a1[0]-a0[0],a1[1]-a0[1])
db = (b1[0]-b0[0],b1[1]-b0[1])
a = (a0[0]+s*da[0],a0[1]+s*da[1])
b = (b0[0]+t*db[0],b0[1]+t*db[1])
d = (b[0]-a[0],b[1]-a[1])
chain = [a,b]
k=(j+1)%m
while k != (i+1)%m:
chain.append(vertices[k]); k=(k+1)%m
left = sum(cross(chain[k],chain[(k+1)%len(chain)]) for k in range(len(chain)))/2
jacobian = -cross(d,da)*cross(d,db)
poly = cross(a,b)*jacobian*(left/area)**(n-2)
total += sum(c/((e[0]+1)*(e[1]+1)) for e,c in poly.dict().items())
return QQ(n*(n-1))/24 * total/area**3
def disk_area(n, precision):
"""Direct signed shoelace-area integral in line-normal coordinates."""
with mp.workdps(precision):
def f(theta):
sn,cs=mp.sin(theta),mp.cos(theta)
fleft=(theta-sn*cs)/mp.pi
return -cs*sn**5*fleft**(n-2)
return + (mp.mpf(4)*n*(n-1)/(3*mp.pi**2)
*mp.quad(f,[0,mp.pi/4,mp.pi/2,3*mp.pi/4,mp.pi]))
def decimal_ball_text(x):
with mp.workdps(180):
return mp.mpf(str(x.center()))
def stored_entries(tid):
doc=numberdb.table(tid)
out={}
# Numbers is a nested parameter dictionary in website documents.
for region, rows in doc['Numbers'].items():
for n, entry in rows.items():
out[(region,int(n))]=entry['number'] if isinstance(entry,dict) else entry
return out
def text_bounds(text):
from decimal import Decimal
from fractions import Fraction
s=str(text)
if '/' in s:
q=Fraction(s);return QQ(q.numerator)/q.denominator, QQ(q.numerator)/q.denominator
if '.' not in s and 'e' not in s.lower():return QQ(s),QQ(s)
d=Decimal(s); q=Fraction(d); ulp=Fraction(10)**d.as_tuple().exponent
return QQ((q-ulp).numerator)/(q-ulp).denominator, QQ((q+ulp).numerator)/(q+ulp).denominator
def main():
if os.environ.get('NUMBERDB_KEY_FROM_STDIN')=='1':
os.environ['NUMBERDB_API_KEY']=sys.stdin.read().strip()
gen=ExpectedHullVolume()
values={(p['region'],int(p['n'])):gen.value(p,100) for p in gen.enumerate()}
bare={k:v['number'] if isinstance(v,dict) else v for k,v in values.items()}
assert not exactness(values)
assert not names_its_rings(os.path.join(os.path.dirname(__file__),'generate.py'))
triangle=[(0,0),(1,0),(0,1)]
square=[(0,0),(1,0),(1,1),(0,1)]
# Controls with answers known before the integration is trusted.
assert polygon_ratio(triangle,3)==QQ(1)/12
assert polygon_ratio(square,3)==QQ(11)/144
assert polygon_ratio(triangle,4)==QQ(1)/6
assert polygon_ratio(square,4)==QQ(11)/72
computed={k:v for k,v in bare.items() if k[0]!='unit-disk'}
independent={k:polygon_ratio(triangle if k[0]=='equilateral-triangle' else square,k[1]) for k in computed}
assert not agrees_with(computed,lambda k:independent[k])
print('Exact polygon edge integrals: 36/36, including normalization controls.',file=sys.stderr)
# An invertible, non-similarity affine map and a translation.
def affine(v):return (2*v[0]+v[1]+3,v[0]+3*v[1]-2)
for k in computed:
vs=triangle if k[0]=='equilateral-triangle' else square
assert polygon_ratio([affine(v) for v in vs],k[1])==independent[k]
print('Exact affine-transformed polygon integrals: 36/36.',file=sys.stderr)
R=RealBallField(600)
pi=R.pi()
special={3:R(35)/(48*pi**2),4:R(35)/(24*pi**2),5:R(175)/(72*pi**2)-R(23023)/(6912*pi**4)}
controls={}
with mp.workdps(180):
# A quadrature control for the same direct edge-area integrator.
ctrl=disk_area(3,170)
assert abs(ctrl-mp.mpf(35)/(48*mp.pi**2))<mp.mpf('1e-160')
max_difference=mp.mpf(0)
for n in range(3,MAX_N+1):
a=disk_area(n,135); b=disk_area(n,175)
assert abs(a-b)<mp.mpf('1e-130')
x=bare[('unit-disk',n)]
discrepancy=abs(decimal_ball_text(x)-b)
assert discrepancy<mp.mpf('1e-118')
max_difference=max(max_difference,discrepancy)
controls[str(n)]={'independent_area_integral':mp.nstr(b,130),'difference':mp.nstr(discrepancy,6)}
if n in special:
assert special[n].is_finite() and R(special[n].rad())<R(10)**-170
assert x.overlaps(special[n])
print('Independent disk signed-area quadratures: 18/18 at 135 and 175 digits.',file=sys.stderr)
for region in REGIONS:
for n in range(3,MAX_N):
assert bare[(region,n)] < bare[(region,n+1)]
for region in REGIONS:
assert (R(bare[(region,4)])-2*R(bare[(region,3)])).contains_zero()
relative_radii=[float(v.rad()/abs(v.center())) for k,v in bare.items() if k[0]=='unit-disk']
report={'entries':len(values),'exact_rationals':len(computed),'certified_disk_enclosures':18,
'measure':measure(values),'max_relative_radius':max(relative_radii),
'polygon_geometric_integrals':36,'affine_checks':36,
'disk_quadrature_precisions':[135,175],'disk_checks':controls,
'known_disk_closed_forms_checked':[3,4,5],
'monotonicity_checks':51,'v4_equals_2v3_checks':3}
if os.environ.get('NUMBERDB_API_KEY'):
stored=stored_entries('T470')
assert set(stored)==set(bare)
with mp.workdps(180):
for k,v in stored.items():
lo,hi=text_bounds(v)
if k in independent: assert lo==hi==independent[k]
else:
b=mp.mpf(controls[str(k[1])]['independent_area_integral'])
assert mp.mpf(str(lo.numerator()))/mp.mpf(str(lo.denominator())) < b < mp.mpf(str(hi.numerator()))/mp.mpf(str(hi.denominator()))
x=bare[k];assert R(lo)<=R(R(x).lower()) and R(R(x).upper())<=R(hi)
# Check the neighboring table's actual stored convex-position probabilities.
neighbor=stored_entries('T468')
for region in REGIONS:
plo,phi=text_bounds(neighbor[(region,4)])
vlo,vhi=text_bounds(stored[(region,3)])
assert max(vlo,(1-phi)/4)<=min(vhi,(1-plo)/4)
report['stored_independent_checks']=54
report['stored_sylvester_relation_checks']=3
for region in REGIONS:
for n in range(3,MAX_N):
assert text_bounds(stored[(region,n)])[1] < text_bounds(stored[(region,n+1)])[0]
report['stored_monotonicity_checks']=51
print('Stored entries and T468 relation checked.',file=sys.stderr)
print(json.dumps(report,indent=2))
if __name__=='__main__':
main()