program.py

back to table · edit · history · where entries came from · files · download

647 bytes, as of the version from 2026-10-29 02:59 (current). Recorded here, not run.

"""A single additional internal simplex angle, using mpmath."""
import mpmath as mp
from fractions import Fraction
mp.mp.dps = 125

def beta(d, k):
    if not (d >= 2 and 0 <= k < d):
        raise ValueError('Require d >= 2 and 0 <= k < d')
    if k == d - 1:
        return Fraction(1, 2)
    if (d, k) == (2, 0):
        return Fraction(1, 6)
    scale = mp.sqrt(2*(d+1))
    cutoff = mp.sqrt(mp.mpf(6)*(d+1)*(mp.mp.dps+30)/(k+1))
    def integrand(x):
        p = (1 + mp.erf(1j*x/scale))/2
        return mp.exp(-x*x/2)*mp.re(p**(d-k))
    return mp.sqrt(2/mp.pi)*mp.quad(integrand, [0, 2, 4, 8, 16, cutoff])

print(mp.nstr(beta(3, 0), 100))