"""tpc_beta 0.1.0 (DRAFT): the Beta layer's arithmetic, shared by the
user score (tpc_sponsor, D160) and the agent score (CR11, another
module).
Reference implementation, Python standard library only. DATA for landing
round 8 of branch tpm_project_closeout; not a declared TPC-PY version and
not TPX code. Every function is a pure function of its arguments and is
written with only +, -, *, /, comparison, sqrt, log and exp, in a fixed
operation order, so that the PHP residence (tpc_beta.php) reproduces every
double bit for bit on one libm and every rendered value byte for byte
through canon() (the tree's %.12g rendering, tpc_strike.py, copied here as
canon_float). Every refused argument raises ValueError with a short fixed
message. Ints are never bools.
Rule sources: CR11 RELIABILITY (D122) "E the mean of Beta(a0 + F, b0 + A)
with the 95 per cent equal-tailed interval and n always shown"; CR12 and
D160.1 (b) "E_u = (a0 + F_u) / (a0 + b0 + F_u + A_u), with the 95 per cent
interval and n_u"; PRIOR a0, b0 = 1, 1 (DECLARED PARAMETERS). The text
fixes the distribution, its mean and the interval's coverage and shape;
the numerical method is this module's (taken by delegation 2026-09-26,
BT-1 to BT-4 below).
beta_mean(a, b) -> float a / (a + b); a, b > 0
beta_interval(a, b, level) -> (lo, hi) equal-tailed: the level
is split half and half; lo = betainc_inv(a, b, (1 - level) / 2),
hi = betainc_inv(a, b, 1 - (1 - level) / 2); 0 < level < 1
betainc(a, b, x) -> float the regularized incomplete
beta function I_x(a, b); a, b > 0; 0 <= x <= 1
betainc_inv(a, b, p) -> float the quantile: the x with
I_x(a, b) = p; 0 <= p <= 1
lgamma(x) -> float log Gamma(x), x > 0
canon_float(x) -> str the tree's %.12g rendering
Method (taken by delegation 2026-09-26; the alternatives in the spec note):
BT-1 lgamma by the Lanczos approximation, g = 7, nine coefficients
(Godfrey's set), reflection below 0.5. Relative error about 1e-15.
BT-2 betainc by the continued fraction of Numerical Recipes (betacf),
evaluated by the modified Lentz method with FPMIN 1e-300, EPS
1e-15 and a fixed iteration cap of 300 (the loop stops early only
on the tolerance and never differs between residences); the
symmetry I_x(a, b) = 1 - I_(1-x)(b, a) used when x > (a + 1) /
(a + b + 2); the prefactor by exp(lgamma(a + b) - lgamma(a) -
lgamma(b) + a log x + b log(1 - x)). Not converged at the cap:
ValueError("betainc").
BT-3 betainc_inv by bisection on [0, 1], exactly 100 halvings, the
midpoint returned; p = 0 gives 0 and p = 1 gives 1.
BT-4 Every value is rendered through canon_float: the two residences
are byte-identical on one host (MEASURED by diff_beta.py); across
hosts a last-ulp libm difference in log or exp can move a twelfth
digit at a rounding boundary (VERIFY on the Windows run).
"""
import math
VERSION = "0.1.0"
LANCZOS_G = 7.0
LANCZOS = (
0.99999999999980993, 676.5203681218851, -1259.1392167224028,
771.32342877765313, -176.61502916214059, 12.507343278686905,
-0.13857109526572012, 9.9843695780195716e-6, 1.5056327351493116e-7)
FPMIN = 1e-300
EPS = 1e-15
MAXIT = 300
INV_ITER = 100
# --------------------------------------------------------------- canon
def canon_float(x):
"""C "%.12g" of a double, as tpc_strike.py's _canon_float builds it:
one rounding to 12 significant digits from "%.11e", fixed notation for
exponents -4..11, exponent notation otherwise, trailing zeros removed;
0 for 0.0 and -0.0."""
if x != x:
raise ValueError("canon")
if x in (float("inf"), float("-inf")):
return "inf" if x > 0 else "-inf"
if x == 0:
return "0"
mant, exp = ("%.11e" % abs(x)).split("e")
digits = mant.replace(".", "")
e = int(exp)
if -4 <= e < 12:
if e >= 0:
ip, fp = digits[:e + 1], digits[e + 1:]
else:
ip, fp = "0", "0" * (-e - 1) + digits
fp = fp.rstrip("0")
out = ip + ("." + fp if fp else "")
else:
rest = digits[1:].rstrip("0")
out = (digits[0] + ("." + rest if rest else "") + "e"
+ ("-" if e < 0 else "+") + "%02d" % abs(e))
return ("-" if x < 0 else "") + out
def _num(x, name):
"""A finite float or int (never a bool), returned as float."""
if type(x) is bool or not isinstance(x, (int, float)):
raise ValueError(name)
x = float(x)
if not math.isfinite(x):
raise ValueError(name)
return x
# -------------------------------------------------------------- lgamma
def lgamma(x):
"""BT-1. log Gamma(x) for x > 0."""
x = _num(x, "lgamma")
if x <= 0.0:
raise ValueError("lgamma")
if x < 0.5:
# reflection: Gamma(x) Gamma(1 - x) = pi / sin(pi x)
return math.log(math.pi / math.sin(math.pi * x)) - lgamma(1.0 - x)
x = x - 1.0
s = LANCZOS[0]
for i in range(1, 9):
s = s + LANCZOS[i] / (x + float(i))
t = x + LANCZOS_G + 0.5
return (0.5 * math.log(2.0 * math.pi) + (x + 0.5) * math.log(t) - t
+ math.log(s))
# ------------------------------------------------------------- betainc
def _betacf(a, b, x):
"""The continued fraction for I_x(a, b) by the modified Lentz method."""
qab = a + b
qap = a + 1.0
qam = a - 1.0
c = 1.0
d = 1.0 - qab * x / qap
if abs(d) < FPMIN:
d = FPMIN
d = 1.0 / d
h = d
m = 1
while m <= MAXIT:
fm = float(m)
m2 = 2.0 * fm
aa = fm * (b - fm) * x / ((qam + m2) * (a + m2))
d = 1.0 + aa * d
if abs(d) < FPMIN:
d = FPMIN
c = 1.0 + aa / c
if abs(c) < FPMIN:
c = FPMIN
d = 1.0 / d
h = h * d * c
aa = -(a + fm) * (qab + fm) * x / ((a + m2) * (qap + m2))
d = 1.0 + aa * d
if abs(d) < FPMIN:
d = FPMIN
c = 1.0 + aa / c
if abs(c) < FPMIN:
c = FPMIN
d = 1.0 / d
de = d * c
h = h * de
if abs(de - 1.0) < EPS:
return h
m += 1
raise ValueError("betainc")
def betainc(a, b, x):
"""BT-2. The regularized incomplete beta function I_x(a, b)."""
a = _num(a, "shape")
b = _num(b, "shape")
x = _num(x, "x")
if a <= 0.0 or b <= 0.0:
raise ValueError("shape")
if x < 0.0 or x > 1.0:
raise ValueError("x")
if x == 0.0:
return 0.0
if x == 1.0:
return 1.0
bt = math.exp(lgamma(a + b) - lgamma(a) - lgamma(b)
+ a * math.log(x) + b * math.log(1.0 - x))
if x < (a + 1.0) / (a + b + 2.0):
return bt * _betacf(a, b, x) / a
return 1.0 - bt * _betacf(b, a, 1.0 - x) / b
def betainc_inv(a, b, p):
"""BT-3. The x with I_x(a, b) = p, by 100 halvings of [0, 1]."""
a = _num(a, "shape")
b = _num(b, "shape")
p = _num(p, "p")
if a <= 0.0 or b <= 0.0:
raise ValueError("shape")
if p < 0.0 or p > 1.0:
raise ValueError("p")
if p == 0.0:
return 0.0
if p == 1.0:
return 1.0
lo = 0.0
hi = 1.0
for _ in range(INV_ITER):
mid = (lo + hi) / 2.0
if betainc(a, b, mid) < p:
lo = mid
else:
hi = mid
return (lo + hi) / 2.0
# ----------------------------------------------------- mean and interval
def beta_mean(a, b):
"""The mean of Beta(a, b): a / (a + b)."""
a = _num(a, "shape")
b = _num(b, "shape")
if a <= 0.0 or b <= 0.0:
raise ValueError("shape")
return a / (a + b)
def beta_interval(a, b, level=0.95):
"""The equal-tailed interval of Beta(a, b) at the level (0.95: the
2.5 and 97.5 per cent quantiles)."""
a = _num(a, "shape")
b = _num(b, "shape")
level = _num(level, "level")
if a <= 0.0 or b <= 0.0:
raise ValueError("shape")
if level <= 0.0 or level >= 1.0:
raise ValueError("level")
alpha = (1.0 - level) / 2.0
return (betainc_inv(a, b, alpha), betainc_inv(a, b, 1.0 - alpha))
# -------------------------------------------------------------- corpus
FNS = ("mean", "interval", "betainc", "betainc_inv", "lgamma", "canon")
def run_case(fn, args):
"""One corpus case: the render, or "ERROR:" and the message."""
if fn not in FNS:
raise ValueError("fn")
try:
if fn == "mean":
return canon_float(beta_mean(args["a"], args["b"]))
if fn == "interval":
lo, hi = beta_interval(args["a"], args["b"], args["level"])
return canon_float(lo) + "|" + canon_float(hi)
if fn == "betainc":
return canon_float(betainc(args["a"], args["b"], args["x"]))
if fn == "betainc_inv":
return canon_float(betainc_inv(args["a"], args["b"], args["p"]))
if fn == "lgamma":
return canon_float(lgamma(args["x"]))
return canon_float(_num(args["x"], "canon"))
except ValueError as ex:
return "ERROR:" + str(ex)
