"""Independent checks for internal simplex angles -- numberdb.org/T598. sage -python checks.py sage -python checks.py --stored # authenticated key needed while a draft The main cross-check uses real Gaussian quadrature for external angles and conic Gauss-Bonnet to recover internal angles. It does not use the generator's imaginary-normal integrand. Closed forms, exact projected Gram matrices and Gram's angle-sum relation supply further geometric checks. Full derivations and the archived skill manifest are recorded in numerical-method.md. """ import json import os import sys import time from functools import lru_cache from fractions import Fraction from math import comb, isqrt from pathlib import Path import numberdb.sage as numberdb from sage.rings.rational_field import QQ from sage.rings.real_arb import RealBallField from sage.rings.complex_arb import ComplexBallField from numberdb._write import to_text from numberdb._generate import _stored_entries from numberdb._compare import _claim from check import exactness, measure, agrees_with, names_its_rings from generate import InternalSimplexAngles, quadrature CHECK_BITS = 630 C = ComplexBallField(CHECK_BITS) R = RealBallField(CHECK_BITS) def sharp(value, digits=100): value = R(value) assert abs(value) < 2, 'Nonfinite or implausible value' assert value.rad_as_ball() < R(10) ** (-digits) * abs(value), 'Imprecise check' return value def overlaps(a, b, digits=100): return sharp(a, digits).overlaps(sharp(b, digits)) @lru_cache(maxsize=None) def external_quadrature(m, q): """g_m(1/q), via a nonnegative real-normal integrand at 630 bits.""" scale = C(2 * q).sqrt() normalization = C(2 * C.pi()).sqrt() tolerance = R(2) ** (-560) cutoff = 34 def f(x, analytic): e = (x / scale).erf() return ((1 + e) ** m + (1 - e) ** m) * (-x * x / 2).exp() / (2 ** m * normalization) z = C.integral(f, 0, cutoff, abs_tol=tolerance, rel_tol=tolerance) assert z.imag().contains_zero() and abs(z.imag()) < tolerance # For real x each Phi power is at most 1; both Gaussian tails. tail = 2 * (-R(cutoff) ** 2 / 2).exp() / (R(cutoff) * (2 * R.pi()).sqrt()) return sharp(z.real().add_error(tail.upper()), 150) @lru_cache(maxsize=None) def external(d, k): if d == k: return QQ(1) if k == 0: return QQ(1) / QQ(d + 1) if k == d - 1: return QQ(1) / QQ(2) return external_quadrature(d - k, k + 1) @lru_cache(maxsize=None) def independent(d, k): """Recover beta using one parity of conic intrinsic volumes.""" if d == k: return R(1) if d == k + 1: return R(QQ(1) / QQ(2)) m = d - k answer = R(QQ(1) / QQ(2)) for j in range(m % 2, m, 2): answer -= comb(m, j) * independent(k + j, k) * external(d, k + j) return sharp(answer, 140) def stored_values(): doc = numberdb.table('T598') written = _stored_entries(doc) values = {} for key, text in written.items(): d, k = map(int, key.split(',')) if '/' in text or ('.' not in text and 'e' not in text.lower()): values[d, k] = QQ(text) else: center, radius = _claim(text) c, r = Fraction(center), Fraction(radius) values[d, k] = R(QQ(c.numerator) / c.denominator).add_error( R(QQ(r.numerator) / r.denominator).upper()) return values, written def run(stored=False): start = time.monotonic() # Known nonzero controls run before any family comparisons. for m, q, known in [(0, 1, QQ(1)), (1, 1, QQ(1)/2), (2, 1, QQ(1)/3), (3, 1, QQ(1)/4)]: assert overlaps(external_quadrature(m, q), known, 145) assert overlaps(quadrature(2, 0), QQ(1)/6, 110) assert overlaps(quadrature(2, 1), QQ(1)/2, 110) print('PASS: 6 nonzero quadrature controls', flush=True) gen = InternalSimplexAngles() keys = [(p['d'], p['k']) for p in gen.enumerate()] assert keys == [(d, k) for d in range(2, 13) for k in range(d)] assert len(keys) == 77 produced = {} written = {} for d, k in keys: value = gen.value({'d':d, 'k':k}, 100) value = value['number'] if isinstance(value, dict) else value produced[d, k] = value written['%d,%d' % (d,k)] = to_text(value, digits=100) assert not exactness(produced) assert not names_its_rings(str(Path(__file__).with_name('generate.py'))) relative_radii = {str(key):float(R(v.rad_as_ball()) / abs(R(v))) for key,v in produced.items() if hasattr(v,'rad_as_ball')} print('PASS: enumeration, exactness and named imports', flush=True) for d, k in keys: high = quadrature(d, k, 140) assert overlaps(produced[d,k], high, 110) print('PASS: all 77 entries agree with 140-digit recomputation', flush=True) for d, k in keys: expected = independent(d,k) assert overlaps(produced[d,k], expected, 110), (d,k) assert not agrees_with(produced, lambda key: independent(*key)) print('PASS: all 77 entries agree with external-angle/Gauss-Bonnet computation', flush=True) values = stored_values()[0] if stored else produced assert list(values) == keys # Written intervals use one unit in their last place, so 97 is the # diagnostic width threshold for values displayed with 100 digits. threshold = 97 if stored else 110 count = {'facets':0, 'codimension_two':0, 'codimension_three':0, 'exact_gram_matrices':0, 'gram_angle_sums':0, 'gauss_bonnet_recurrences':0} for d,k in keys: v = sharp(values[d,k], threshold) assert v > 0 and v < 1 assert overlaps(v, independent(d,k), threshold) m=d-k # Project e_j - e_0 orthogonally off the face direction space: # u_j = e_j - mean(e_0,...,e_k). Build their actual coordinates. vectors = [[(QQ(1) if h==j else QQ(0)) - (QQ(1)/QQ(k+1) if h<=k else QQ(0)) for h in range(d+1)] for j in range(k+1,d+1)] gram = [[sum(x*y for x,y in zip(u,w)) for w in vectors] for u in vectors] inv = [[QQ(int(i==j))-QQ(1)/QQ(d+1) for j in range(m)] for i in range(m)] for i in range(m): for j in range(m): assert gram[i][j] == QQ(int(i==j))+QQ(1)/QQ(k+1) assert sum(gram[i][h]*inv[h][j] for h in range(m)) == int(i==j) count['exact_gram_matrices']+=1 if m==1: assert overlaps(v, QQ(1)/2, threshold); count['facets']+=1 if m==2: assert overlaps(v, (R(1)/d).arccos()/(2*R.pi()), threshold) count['codimension_two']+=1 if m==3: assert overlaps(v, QQ(1)/8 - 3*(R(1)/d).arcsin()/(4*R.pi()), threshold) count['codimension_three']+=1 # Evaluate the full parity identity using the chosen values themselves. total=R(v) for j in range(m%2,m,2): inner=R(1) if j==0 else (R(QQ(1)/2) if j==1 else R(values[k+j,k])) total += comb(m,j)*inner*external(d,k+j) assert (total-QQ(1)/2).contains_zero() and total.rad_as_ball()