version_h.py

1import mpmath as mp2import numpy as np3import sympy as sp45TOP = 76MC_DIMS = (3, 4, 5)7MC_SAMPLES = 10**88BATCH = 10**6910def beta_sym(x, y):11    return sp.gamma(x) * sp.gamma(y) / sp.gamma(x + y)1213def exact_miss(d):14    th, ps = sp.symbols("th ps", positive=True)15    alpha = sp.Rational(d * d - 1, 2)16    inner = sp.integrate(sp.sin(th) ** (d * d), (th, 0, sp.pi / 2 - ps))17    body = sp.expand_trig(sp.expand(inner * sp.sin(ps) ** (d - 2)))18    num = 2 * sp.integrate(body, (ps, 0, sp.pi / 2))19    denom = beta_sym(sp.Rational(1, 2), sp.Rational(d - 1, 2)) * beta_sym(20        sp.Rational(1, 2), alpha + 121    )22    return sp.radsimp(sp.expand(sp.nsimplify(sp.simplify(2**d * num / denom))))2324def quad_miss(d, dps):25    saved = mp.mp.dps26    mp.mp.dps = dps27    alpha = mp.mpf(d * d - 1) / 228    half = mp.mpf(1) / 22930    def shell(h):31        return mp.betainc(mp.mpf(d - 1) / 2, half, 0, h * h, regularized=True)3233    num = mp.quad(lambda h: (1 - h * h) ** alpha * shell(h), [0, half, 1])34    out = mp.mpf(2) ** d * num / mp.beta(half, alpha + 1)35    mp.mp.dps = saved36    return out3738def half_ball(rng, count, dim):39    pts = rng.standard_normal((count, dim, dim))40    pts /= np.linalg.norm(pts, axis=2, keepdims=True)41    pts *= rng.random((count, dim, 1)) ** (1.0 / dim)42    pts[:, :, dim - 1] = np.abs(pts[:, :, dim - 1])43    return pts4445def flat_face_rate(rng, dim, samples):46    hits = 047    done = 048    while done < samples:49        take = min(BATCH, samples - done)50        pts = half_ball(rng, take, dim)51        mat = np.concatenate(52            [pts[:, :, : dim - 1], np.ones((take, dim, 1))], axis=253        )54        rhs = pts[:, :, dim - 1][:, :, None]55        sol = np.linalg.solve(mat, rhs)[:, :, 0]56        slope = np.linalg.norm(sol[:, : dim - 1], axis=1)57        hits += int(np.count_nonzero(np.abs(sol[:, dim - 1]) <= slope))58        done += take59    return hits / samples6061def main():62    mp.mp.dps = 3063    print("VERSION H: d uniform points, the hyperplane through them misses the flat face")64    print("  reduction: unoriented normal n with n_d = t >= 0 and signed offset h")65    print("  the section is a (d-1)-ball of radius sqrt(1-h^2) and stays above the")66    print("  flat face exactly when h >= sqrt(1-t^2)")67    print("  the within-hyperplane simplex integral scales as radius^(d^2-1)")68    print("  the half-ball Cartesian d-tuple measure is 2^-d of the full-ball one")6970    forms = {}71    print("  exact values by symbolic reduction")72    for d in range(2, TOP + 1):73        p = exact_miss(d)74        forms[d] = p75        print(f"  d = {d}  {p}")7677    print("  decimals to ten places")78    for d in range(2, TOP + 1):79        print(f"  d = {d}  {float(sp.N(forms[d], 30)):.10f}")8081    print("  quadrature against the exact values, 60 and 80 working digits")82    for d in range(2, TOP + 1):83        lo = quad_miss(d, 60)84        hi = quad_miss(d, 80)85        mp.mp.dps = 8086        target = mp.mpf(str(sp.N(forms[d], 90)))87        print(88            f"  d = {d}  60 against 80 digits {mp.nstr(abs(lo - hi), 3)}"89            f"  quadrature against exact {mp.nstr(abs(hi - target), 3)}"90        )91        mp.mp.dps = 309293    print("  parity read off the exact values")94    for d in range(2, TOP + 1):95        powers = sorted(96            int(term.as_coeff_exponent(sp.pi)[1]) for term in sp.Add.make_args(forms[d])97        )98        shape = " + ".join("Q" if e == 0 else f"Q Pi^{e}" for e in powers)99        print(f"  d = {d}  pi powers {powers}  {shape}")100101    rng = np.random.default_rng(20250828)102    print(f"  random-point check, {MC_SAMPLES} samples per dimension")103    for d in MC_DIMS:104        flat_exact = 1 - mp.mpf(str(sp.N(forms[d], 30)))105        est = mp.mpf(flat_face_rate(rng, d, MC_SAMPLES))106        se = mp.sqrt(est * (1 - est) / MC_SAMPLES)107        print(108            f"  d = {d}  flat-face estimate {mp.nstr(est, 8)}"109            f"  exact {mp.nstr(flat_exact, 8)}"110            f"  deviation {mp.nstr(est - flat_exact, 3)}  one sigma {mp.nstr(se, 3)}"111        )112113main()