The Principia Machina/TPC-PY/tpc beta.py: Difference between revisions

Latest revision as of 04:37, 28 September 2026

"""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)