"""Independent checks of expected hull-vertex counts -- numberdb.org/T471. Run sage -python checks.py for local tests; set NUMBERDB_CHECK_STORED=1 and NUMBERDB_API_KEY to check the private draft and its two companion tables. NUMBERDB_KEY_FROM_STDIN=1 also accepts the key from stdin. For a convex polygon of area A, an unordered pair of sample points forms a hull edge iff the other n-2 points lie on one side of its supporting line. Changing the two points to line coordinates gives L^3/3 after integrating along a chord of length L. Parameterizing the chord by its intersections P(s), Q(t) with two boundary edges cancels L^3 against the line-coordinate Jacobian. The remaining weight is abs(cross(Q-P,P')*cross(Q-P,Q')). Multiply its integral of (F^(n-2)+(1-F)^(n-2)) by binom(n,2)/(3*A^2). In the right unit triangle the three edge pairs have weight uv and cap fraction uv. In the unit square the four adjacent pairs have weight uv and cap fraction uv/2, and the two opposite pairs have weight 1 and cap fraction (u+v)/2. The routines below integrate their binomial expansions exactly; they do not use harmonic sums or the generator's code. For the disk an independent area integral is taken from Efron's technical report SOL-ONR-103, printed p.15. The conditional chord variance is (1-p^2)/3 and the projection density is 2*sqrt(1-p^2)/pi. Thus v_D(m) equals 8*binom(m,2)/(3*pi^2) times the integral from -1 to 1 of p*(1-p^2)^2*G(p)^(m-2), where G(p)=(acos(-p)+p*sqrt(1-p^2))/pi. Efron's identity then gives e_D(n)=n*(1-v_D(n-1)). This differs from the expected-edge integral used by the generator. mpmath quadrature at 130 and 170 decimal digits is an agreement check, not the source of certification. The known three-point area 35/(48*pi^2) is tested before other quadratures. All output can be redirected to a log; no credentials are printed. """ import os import sys import json from math import comb from decimal import Decimal from fractions import Fraction from pathlib import Path import numberdb.sage as numberdb from sage.rings.rational_field import QQ from sage.rings.real_arb import RealBallField from generate import ExpectedHullVertices, REGIONS, ellipse from check import exactness, measure, agrees_with, names_its_rings from numberdb._write import to_text import mpmath as mp def triangle_edges(n): m = n - 2 integral = QQ(1) / n**2 integral += sum((QQ((-1)**k * comb(m,k)) / (k+2)**2 for k in range(m+1)), QQ(0)) return 4 * comb(n,2) * integral def square_edges(n): m = n - 2 adjacent = QQ(1) / (2**m * (m+2)**2) adjacent += sum((QQ((-1)**k * comb(m,k)) / (2**k * (k+2)**2) for k in range(m+1)), QQ(0)) opposite = QQ(2*(2**(m+2)-2)) / (2**m*(m+1)*(m+2)) return QQ(comb(n,2))/3 * (4*adjacent + 2*opposite) def disk_area(m, dps): with mp.workdps(dps): def integrand(p): G = (mp.acos(-p)+p*mp.sqrt(1-p*p))/mp.pi return p*(1-p*p)**2 * G**(m-2) return 8*mp.mpf(comb(m,2))/(3*mp.pi**2)*mp.quad(integrand, [-1,0,1]) def disk_edges_from_area(n,dps): with mp.workdps(dps): return n*(1-disk_area(n-1,dps)) def rational_decimal(text): q = Fraction(Decimal(text)) return QQ(q.numerator)/q.denominator def check_local(): # Controls come first, including deliberately distinct known expectations. assert triangle_edges(3) == square_edges(3) == 3 assert triangle_edges(4) == QQ(11)/3 assert square_edges(4) == QQ(133)/36 R = RealBallField(600) control = ellipse(3) assert control.is_finite() and R(control.rad()) < R(10)**-115 assert (control-3).contains_zero() with mp.workdps(170): for prec in [130,170]: area = disk_area(3,prec) assert mp.isfinite(area) assert abs(area-35/(48*mp.pi**2)) < mp.mpf('1e-125') print('PASS known controls: polygon n=3,4; disk n=3; disk triangle-area integral',flush=True) g = ExpectedHullVertices() values = {(p['region'],int(p['n'])):g.value(p,100) for p in g.enumerate()} raw = {k:(v['number'] if isinstance(v,dict) else v) for k,v in values.items()} assert not exactness(values) assert not names_its_rings(str(Path(__file__).with_name('generate.py'))) polygons = {k:v for k,v in raw.items() if k[0]!='unit-disk'} mismatch = agrees_with(polygons,lambda k:triangle_edges(k[1]) if k[0]=='equilateral-triangle' else square_edges(k[1])) assert not mismatch, mismatch print('PASS all 34 rationals against exact hull-edge integrals',flush=True) quad = {} max_rad = R(0) for n in range(4,21): with mp.workdps(180): a = disk_edges_from_area(n,130) b = disk_edges_from_area(n,170) assert mp.isfinite(a) and mp.isfinite(b) assert abs(a-b) < mp.mpf('1e-125') v = raw['unit-disk',n] assert v.is_finite() and R(v.rad()) < R(10)**-115 # mpmath is a diagnostic comparison, with 115-digit tolerance. mid = mp.mpf(str(v.center())) assert abs(mid-b) < mp.mpf('1e-115') quad[n] = mp.nstr(b,175) max_rad = max(max_rad,R(v.rad())) print('PASS all 17 disk values against independent signed-area integrals at 130 and 170 digits',flush=True) assert (R(raw['unit-disk',4]) - (4-R(35)/(12*R.pi()**2))).contains_zero() for region in REGIONS: for n in range(4,21): v=raw[region,n] assert v > 3 and v < n print('PASS bounds 3 < e_K(n) < n and disk four-point closed form',flush=True) print('MEASURE',json.dumps(measure(values),default=str),flush=True) print('MAX_DISK_RADIUS',str(max_rad),flush=True) print('NUMERICAL_EVIDENCE',json.dumps([ {'region':k[0],'n':k[1],'number':to_text(v,100), 'radius':str(v.rad()) if k[0]=='unit-disk' else '0', 'independent':quad.get(k[1]) if k[0]=='unit-disk' else str(triangle_edges(k[1]) if k[0]=='equilateral-triangle' else square_edges(k[1]))} for k,v in raw.items()]),flush=True) return raw,quad def flatten(doc): return {(region,int(n)):entry['number'] if isinstance(entry,dict) else entry for region,rows in doc['Numbers'].items() for n,entry in rows.items()} def text_interval(s): s=str(s) if '/' in s or ('.' not in s and 'e' not in s.lower()): q=QQ(s) return q,q q=rational_decimal(s) exp=Decimal(s).as_tuple().exponent unit=QQ(10)**exp return q-unit,q+unit def check_stored(raw,quad): g=ExpectedHullVertices() report=g.verify(sample=None) print('VERIFY',report,flush=True) assert report.ok # check.stored requires Django; the public client supplies the same entries. stored=flatten(numberdb.table('T471')) volumes=flatten(numberdb.table('T470')) probs=flatten(numberdb.table('T468')) assert set(stored)==set(raw) for (region,n),text in stored.items(): lo,hi=text_interval(text) if region!='unit-disk': expected=triangle_edges(n) if region=='equilateral-triangle' else square_edges(n) assert lo==hi==expected else: assert lo < rational_decimal(quad[n]) < hi assert lo>3 and hi