version_l.py

1from fractions import Fraction23import mpmath as mp4import numpy as np56TOP = 117MC_TOP = 78MC_SAMPLES = 10**79BATCH = 10**61011def wallis_table(top):12    table = {}13    even = Fraction(16, 3)14    odd = Fraction(3, 8)15    table[2] = ("r/Pi^2", even)16    table[3] = ("rational", odd)17    k = 118    while 2 * k + 2 <= top or 2 * k + 3 <= top:19        even *= Fraction(4 * k * (k + 1), (2 * k + 1) * (2 * k + 3))20        odd *= Fraction((2 * k + 1) * (2 * k + 3), (2 * k + 2) * (2 * k + 4))21        if 2 * k + 2 <= top:22            table[2 * k + 2] = ("r/Pi^2", even)23        if 2 * k + 3 <= top:24            table[2 * k + 3] = ("rational", odd)25        k += 126    return table2728def value(form, rational):29    if form == "r/Pi^2":30        return mp.mpf(rational.numerator) / mp.mpf(rational.denominator) / mp.pi**231    return mp.mpf(rational.numerator) / mp.mpf(rational.denominator)3233def half_ball(rng, count, dim):34    pts = rng.standard_normal((count, dim))35    pts /= np.linalg.norm(pts, axis=1, keepdims=True)36    pts *= rng.random((count, 1)) ** (1.0 / dim)37    pts[:, dim - 1] = np.abs(pts[:, dim - 1])38    return pts3940def flat_face_rate(rng, dim, samples):41    hits = 042    done = 043    while done < samples:44        take = min(BATCH, samples - done)45        p = half_ball(rng, take, dim)46        q = half_ball(rng, take, dim)47        gap = p[:, dim - 1] - q[:, dim - 1]48        good = gap != 0.049        step = np.where(good, p[:, dim - 1] / np.where(good, gap, 1.0), 0.0)50        cross = p[:, : dim - 1] + step[:, None] * (q[:, : dim - 1] - p[:, : dim - 1])51        hits += int(np.count_nonzero(good & (np.linalg.norm(cross, axis=1) <= 1.0)))52        done += take53    return hits / samples5455def main():56    mp.mp.dps = 3057    table = wallis_table(TOP)58    print("VERSION L: two uniform points, the line through them meets the flat face")59    print("  seeds f(2) = 16/(3 Pi^2) and f(3) = 3/8")60    print("  even step f(2k+2)/f(2k) = 4k(k+1)/((2k+1)(2k+3))")61    print("  odd step  f(2k+3)/f(2k+1) = (2k+1)(2k+3)/((2k+2)(2k+4))")62    print(f"  exact rationals to d = {TOP}")63    for d in sorted(table):64        form, r = table[d]65        v = value(form, r)66        shown = f"({r})/Pi^2" if form == "r/Pi^2" else f"{r}"67        print(f"  d = {d:>2}  {shown:>26}  {form:>9}  {mp.nstr(v, 12)}")6869    print("  page decimals")70    for d, digits in ((2, 7), (3, 3), (4, 7), (5, 6), (6, 7), (7, 10)):71        form, r = table[d]72        print(f"  d = {d}  {mp.nstr(value(form, r), digits)}")7374    rng = np.random.default_rng(20250828)75    print(f"  random-point check, {MC_SAMPLES} samples per dimension")76    for d in range(2, MC_TOP + 1):77        form, r = table[d]78        exact = value(form, r)79        est = mp.mpf(flat_face_rate(rng, d, MC_SAMPLES))80        se = mp.sqrt(est * (1 - est) / MC_SAMPLES)81        print(82            f"  d = {d}  estimate {mp.nstr(est, 8)}  exact {mp.nstr(exact, 8)}"83            f"  deviation {mp.nstr(est - exact, 3)}  one sigma {mp.nstr(se, 3)}"84        )8586main()