Follow the source behind the result.

The registered Python source is shown as text. Displaying it does not execute the code.

"""v100 -- the numerology null test: look-elsewhere-corrected hit statistics.

Executes the long-standing TODO from the legacy CMB paper ("Numerology grammar
test: enumerate every reasonable algebraic combination of the primitives;
count how many also land within the observed windows") as a machine check of
the active suite.  It answers the only statistically meaningful form of the
question "is TFPT numerology?":

    Under an EXPLICIT null model -- random formulas drawn from a declared
    grammar that is at least as expressive as TFPT's own formula shapes,
    complexity-matched per observable -- what is the probability that a
    random "theory" reproduces the TFPT data scorecard?

Four layers:

  (A) EXACT LOOK-ELSEWHERE CENSUS.  For each scored observable the full
      grammar class is counted (no sampling): hit density p_i = (#expressions
      inside the observable's honest data window) / (#expressions inside the
      declared two-decade plausibility window).  The joint formula-fishing
      probability is prod p_i.  Robustness: the census is repeated with the
      coefficient bound and the plausibility window varied; the answer must
      be stable in order of magnitude.
  (B) MONTE CARLO over random dial settings (deterministic seeds): random
      rational coefficients inside TFPT's own formula skeletons -- with the
      seed phi0 fixed (MC-A) and with a random seed phi* drawn log-uniformly
      per pseudo-theory (MC-B).  Reported: hit-count distribution, maximum
      hits over all trials (TFPT scores 13/13 by construction).
  (C) NEGATIVE CONTROLS (the test must have power): perturbing the seed phi0
      by 1%/10%/30% collapses TFPT's own hit count; shuffling which measured
      value is assigned to which formula collapses it to ~0.
  (D) THE ALPHA CUBIC CLASS: all complexity-matched variants of the F_U1
      fixed-point equation (coefficients, M, transport exponent, resolvent
      power) are root-solved; only the TFPT equation lands within the
      achieved 4e-8 window of CODATA alpha^-1 -> p_alpha <= hits/N_class.

HONEST SCOPE (do not oversell):
  * No statistical test can give absolute certainty; the result is a
    look-elsewhere-corrected null probability CONDITIONAL on the declared
    grammar.  The grammar is published in full in the CONFIG block below and
    provably contains every TFPT formula it scores (fairness checks).
  * Data windows are conservative: each window is the MAXIMUM of TFPT's
    achieved deviation, the documented experimental sigma, the quote
    resolution of the documented central value, and (for quark ratios) the
    5% scheme-spread tolerance of v18.  Tension observables (theta13, s23,
    delta_CKM, m_mu/m_tau, beta) therefore enter with their honest LARGE
    windows and contribute almost nothing to the significance.
  * Measured centrals are exactly the ones already documented in this repo
    (predictions_frozen.json "experiment" fields, v62, v88, v18, tfpt_2
    tables); nothing is imported fresh for this test.
  * This is a numerical/statistical module (floats, RNG, numpy); it is
    Python-only by the suite's convention (like v62/v64/v65) and flagged as
    such in the Wolfram README.  Complementary to v84 (pre-registration
    freeze): v84 prevents future look-elsewhere drift, v100 quantifies the
    retrospective look-elsewhere burden.
"""
import json
import math
import os
import random
from fractions import Fraction

import numpy as np
import mpmath as mp

from tfpt_constants import check, summary, reset, phi0 as phi0_mp, c3 as c3_mp

# ---------------------------------------------------------------- CONFIG ----
PHI0 = float(phi0_mp)                    # retained seed (from c3 = 1/(8 pi))
LAM = math.sqrt(PHI0 * (1.0 - PHI0))     # Cabibbo lambda_C(phi0)
PI = math.pi
C3 = float(c3_mp)

# one-term grammar G1: value = (p/q) * phi^a * lam^c * (1+lam)^d * pi^b * e^(u/w)
G1_PQ_MAX = 130          # covers TFPT's largest coefficient 55/117
G1_A = range(-2, 3)      # phi exponent
G1_C = range(-3, 4)      # lambda exponent (s13 needs lam^3)
G1_D = range(-1, 2)      # (1+lam) exponent (s23 needs -1)
G1_B = range(-1, 2)      # pi exponent (beta needs -1)
G1_EXP_NUM, G1_EXP_DEN = 6, 6   # e^(u/w), |u|<=6, w<=6 (theta13 needs -5/6)

# two-term grammar G2: value = (p1/q1)*B1 +/- (p2/q2)*B2,
#   B = phi^a * lam^c * pi^b with a,c in 0..2, b in -1..1, coefficients <= 12
G2_PQ_MAX = 12
G2_A = range(0, 3)
G2_C = range(0, 3)
G2_B = range(-1, 2)

DECADE = 10.0            # plausibility window [m/DECADE, m*DECADE]
MC_TRIALS = 200_000
RNG_SEED = 20260611


def _half_ulp(s):
    """Half a unit-in-last-place of a documented decimal quote string."""
    if "." not in s:
        return 0.5
    return 0.5 * 10.0 ** (-len(s.split(".")[1]))


# ------------------------------------------------------- scored observables ----
# Measured centrals/sigmas are the repo-documented ones (see docstring).
# fields: id, predict(phi,lam), measured-quote, sigma (or None), extra window
# floor (quark scheme spread, v18), arity (1 or 2), TFPT G1/G2 encoding.
def _obs_table():
    return [
        # ---- two-term shapes (G2) ----
        dict(id="SIN2_THETA12", arity=2, quote="0.307", sigma=0.012, floor=0.0,
             pred=lambda f, l: 1.0 / 3.0 - f / 2.0,
             g2=((1, 3, 0, 0, 0), -1, (1, 2, 1, 0, 0)),
             src="NuFIT 6.0 (v62)"),
        dict(id="OMEGA_B", arity=2, quote="0.0493", sigma=None, floor=0.0,
             pred=lambda f, l: f - (f / 4.0) / PI,
             g2=((1, 1, 1, 0, 0), -1, (1, 4, 1, 0, -1)),
             src="Planck 2018 (registry)"),
        dict(id="DELTA_CKM_RAD", arity=2, quote="1.146681", sigma=3.0 * PI / 180.0,
             floor=0.0, cap=2.0 * PI,
             pred=lambda f, l: PI / 3.0 + 3.0 * l * l,
             g2=((1, 3, 0, 0, 1), +1, (3, 1, 0, 2, 0)),
             src="gamma_PDG 65.7 +- 3.0 deg (v88)"),
        # ---- one-term shapes (G1) ----
        dict(id="SIN2_THETA13", arity=1, quote="0.02195", sigma=0.00058, floor=0.0,
             pred=lambda f, l: f * math.exp(-5.0 / 6.0),
             g1=(1, 1, 1, 0, 0, 0, -5, 6), src="NuFIT 6.0 (v62)"),
        dict(id="BETA_RAD", arity=1, quote="0.00375246", sigma=0.074 * PI / 180.0,
             floor=0.0,
             pred=lambda f, l: f / (4.0 * PI),
             g1=(1, 4, 1, 0, 0, -1, 0, 1), src="ACT DR6 0.215+-0.074 deg (v62)"),
        dict(id="LAMBDA_C", arity=1, quote="0.2245", sigma=None, floor=0.0,
             pred=lambda f, l: l,
             g1=(1, 1, 0, 1, 0, 0, 0, 1), src="PDG (tfpt_2 CKM table)"),
        dict(id="S23_CKM", arity=1, quote="0.041", sigma=None, floor=0.0,
             pred=lambda f, l: f / (1.0 + l),
             g1=(1, 1, 1, 0, -1, 0, 0, 1), src="PDG ~0.041 (registry)"),
        dict(id="S13_CKM", arity=1, quote="0.0038", sigma=None, floor=0.0,
             pred=lambda f, l: l ** 3 / 3.0,
             g1=(1, 3, 0, 3, 0, 0, 0, 1), src="PDG ~0.0038 (registry)"),
        dict(id="MMU_OVER_MTAU", arity=1, quote="0.05946", sigma=None, floor=0.0,
             pred=lambda f, l: 8.0 / 7.0 * f,
             g1=(8, 7, 1, 0, 0, 0, 0, 1), src="PDG (registry)"),
        dict(id="ME_OVER_MMU", arity=1, quote="0.004836", sigma=None, floor=0.0,
             pred=lambda f, l: 12.0 / 7.0 * f * f,
             g1=(12, 7, 2, 0, 0, 0, 0, 1), src="PDG (registry)"),
        dict(id="MU_OVER_MD", arity=1, quote="0.47", sigma=None, floor=0.05,
             pred=lambda f, l: 55.0 / 117.0,
             g1=(55, 117, 0, 0, 0, 0, 0, 1), src="lattice/FLAG ~0.47 (v18 tol 5%)"),
        dict(id="MC_OVER_MS", arity=1, quote="13.6", sigma=None, floor=0.05,
             pred=lambda f, l: 34.0 / 47.0 / f,
             g1=(34, 47, -1, 0, 0, 0, 0, 1), src="PDG MSbar ~13.6 (v18 tol 5%)"),
        dict(id="MT_OVER_MB", arity=1, quote="41", sigma=None, floor=0.05,
             pred=lambda f, l: 3.0 / 26.0 / (f * f),
             g1=(3, 26, -2, 0, 0, 0, 0, 1), src="PDG ~41 (registry)"),
    ]


def build_observables():
    """Attach measured central m and conservative relative window w to each."""
    obs = _obs_table()
    for o in obs:
        m = float(o["quote"])
        pred = o["pred"](PHI0, LAM)
        achieved = abs(pred - m) / m
        sig = (o["sigma"] / m) if o["sigma"] else 0.0
        res = _half_ulp(o["quote"]) / m
        o["m"] = m
        o["pred_val"] = pred
        o["window"] = max(achieved, sig, res, o["floor"])
    return obs


# ------------------------------------------------ layer A: exact census ----
def g1_structures(phi, lam):
    """All G1 structure factors F = phi^a lam^c (1+lam)^d pi^b e^(u/w)."""
    exps = [(0, 1)] + [(u, w) for w in range(1, G1_EXP_DEN + 1)
                       for u in range(-G1_EXP_NUM, G1_EXP_NUM + 1)
                       if u != 0 and math.gcd(abs(u), w) == 1]
    fs = []
    for a in G1_A:
        for c in G1_C:
            for d in G1_D:
                for b in G1_B:
                    base = phi ** a * lam ** c * (1 + lam) ** d * PI ** b
                    for (u, w) in exps:
                        fs.append(base * math.exp(u / w))
    return np.array(fs)


def census_one_term(F, m, window, decade, pq_max, cap=None):
    """Count G1 expressions (p/q)*F inside the hit / plausibility windows."""
    lo_h, hi_h = m * (1 - window), m * (1 + window)
    lo_d, hi_d = m / decade, m * decade
    if cap is not None:
        hi_d = min(hi_d, cap)
    hits = dec = 0
    for q in range(1, pq_max + 1):
        for lo, hi, acc in ((lo_h, hi_h, "h"), (lo_d, hi_d, "d")):
            plo = np.ceil(q * lo / F)
            phi_ = np.floor(q * hi / F)
            cnt = np.clip(np.minimum(phi_, pq_max) - np.maximum(plo, 1) + 1,
                          0, None).sum()
            if acc == "h":
                hits += cnt
            else:
                dec += cnt
    return float(hits), float(dec)


def g2_terms(phi, lam, pq_max):
    """All G2 term values r*B (sorted numpy array)."""
    rats = [p / q for q in range(1, pq_max + 1) for p in range(1, pq_max + 1)
            if math.gcd(p, q) == 1]
    vals = []
    for a in G2_A:
        for c in G2_C:
            for b in G2_B:
                B = phi ** a * lam ** c * PI ** b
                vals.extend(r * B for r in rats)
    return np.sort(np.array(vals))


def census_two_term(T, m, window, decade, cap=None):
    """Count G2 expressions t1 +/- t2 inside the hit / plausibility windows."""
    lo_h, hi_h = m * (1 - window), m * (1 + window)
    lo_d, hi_d = m / decade, m * decade
    if cap is not None:
        hi_d = min(hi_d, cap)

    def count(lo, hi):
        # plus branch: t2 in [lo - t1, hi - t1]
        n = (np.searchsorted(T, hi - T, side="right")
             - np.searchsorted(T, lo - T, side="left")).sum()
        # minus branch: t2 in [t1 - hi, t1 - lo]
        n += (np.searchsorted(T, T - lo, side="right")
              - np.searchsorted(T, T - hi, side="left")).sum()
        return float(n)

    return count(lo_h, hi_h), count(lo_d, hi_d)


def run_census(obs, pq1, pq2, decade):
    F = g1_structures(PHI0, LAM)
    T = g2_terms(PHI0, LAM, pq2)
    out = {}
    for o in obs:
        cap = o.get("cap")
        if o["arity"] == 1:
            h, d = census_one_term(F, o["m"], o["window"], decade, pq1, cap)
        else:
            h, d = census_two_term(T, o["m"], o["window"], decade, cap)
        out[o["id"]] = (h, d, h / d if d else 1.0)
    return out


# ----------------------------------------- layer B: Monte Carlo pseudo-theories ----
def _rat(rng, pmax):
    return (int(rng.random() * pmax) + 1) / (int(rng.random() * pmax) + 1)


def random_theory_values(rng, phi, lam):
    """Random rational dials inside TFPT's own 13 formula skeletons."""
    eu = int(rng.random() * (2 * G1_EXP_NUM + 1)) - G1_EXP_NUM
    ew = int(rng.random() * G1_EXP_DEN) + 1
    return [
        _rat(rng, 12) - _rat(rng, 12) * phi,                  # theta12
        _rat(rng, 12) * phi - _rat(rng, 12) * phi / PI,       # Omega_b
        _rat(rng, 12) * PI + _rat(rng, 12) * lam * lam,       # delta_CKM
        _rat(rng, 12) * phi * math.exp(eu / ew),              # theta13
        _rat(rng, 12) * phi / PI,                             # beta (rad)
        _rat(rng, 12) * lam,                                  # lambda_C
        _rat(rng, 12) * phi / (1.0 + lam),                    # s23
        _rat(rng, 12) * lam ** 3,                             # s13
        _rat(rng, 12) * phi,                                  # mmu/mtau
        _rat(rng, 12) * phi * phi,                            # me/mmu
        _rat(rng, 130),                                       # mu/md
        _rat(rng, 130) / phi,                                 # mc/ms
        _rat(rng, 130) / (phi * phi),                         # mt/mb
    ]


MC_ORDER = ["SIN2_THETA12", "OMEGA_B", "DELTA_CKM_RAD", "SIN2_THETA13",
            "BETA_RAD", "LAMBDA_C", "S23_CKM", "S13_CKM", "MMU_OVER_MTAU",
            "ME_OVER_MMU", "MU_OVER_MD", "MC_OVER_MS", "MT_OVER_MB"]


def hit_count(values, obs_by_id):
    k = 0
    for oid, v in zip(MC_ORDER, values):
        o = obs_by_id[oid]
        if abs(v - o["m"]) <= o["window"] * o["m"]:
            k += 1
    return k


def run_mc(obs_by_id, trials, seed, random_seed_constant):
    rng = random.Random(seed)
    counts = [0] * 14
    for _ in range(trials):
        if random_seed_constant:
            phi = math.exp(rng.random() * (math.log(0.5) - math.log(0.005))
                           + math.log(0.005))
        else:
            phi = PHI0
        lam = math.sqrt(phi * (1.0 - phi))
        counts[hit_count(random_theory_values(rng, phi, lam), obs_by_id)] += 1
    return counts


# --------------------------------------------- layer D: the alpha cubic class ----
def alpha_class_census():
    """Root-solve every complexity-matched F_U1 variant; count CODATA hits.

    TFPT instance: r1=2, r2=4/5, M=41, transport e^(-2a), resolvent power -5/4.
    """
    pb = 1.0 / (6.0 * PI)
    dt = 48.0 * C3 ** 4
    r1s = np.array([1.0, 2.0, 3.0, 4.0, 5.0])
    r2s = np.array(sorted({float(Fraction(p, q))
                           for p in range(1, 11) for q in range(1, 6)}))
    Ms = np.arange(1, 61, dtype=float)
    ts = np.array([1.0, 2.0, 3.0])
    pws = np.array([-1.25, -1.0, -0.5])

    R1, R2, M, T, PW = np.meshgrid(r1s, r2s, Ms, ts, pws, indexing="ij")
    R1, R2, M, T, PW = (x.ravel() for x in (R1, R2, M, T, PW))
    n = R1.size

    lo = np.full(n, 1e-4)
    hi = np.full(n, 0.2)

    def f(a):
        Q = dt * np.exp(-T * a)
        ps = pb + Q * (1.0 - Q) ** PW
        return a ** 3 - R1 * C3 ** 3 * a ** 2 - R2 * C3 ** 6 * M * np.log(1.0 / ps)

    flo = f(lo)
    for _ in range(70):
        mid = 0.5 * (lo + hi)
        fm = f(mid)
        left = (flo * fm) <= 0
        hi = np.where(left, mid, hi)
        lo = np.where(left, lo, mid)
        flo = np.where(left, flo, fm)
    roots = 0.5 * (lo + hi)
    ainv = 1.0 / roots

    codata = 137.035999177          # CODATA 2022 (v62)
    window = 4.0e-8                 # TFPT's achieved |Delta(alpha^-1)|, >= sigma
    hits = int(np.sum(np.abs(ainv - codata) <= window))

    # fairness: the TFPT combination is in the grid and reproduces the root
    sel = ((R1 == 2.0) & (np.abs(R2 - 0.8) < 1e-12) & (M == 41.0)
           & (T == 2.0) & (PW == -1.25))
    tfpt_ainv = float(ainv[sel][0])
    return n, hits, tfpt_ainv


# ----------------------------------------------------------------- driver ----
def run():
    reset()
    print("v100 numerology null test (look-elsewhere census + MC + controls)")

    obs = build_observables()
    obs_by_id = {o["id"]: o for o in obs}

    # (0) registry lock: the scored predictions are the frozen ones
    reg_path = os.path.join(os.path.dirname(os.path.abspath(__file__)),
                            "predictions_frozen.json")
    with open(reg_path) as fh:
        reg = {e["id"]: e["frozen_value"] for e in json.load(fh)["predictions"]}
    reg_map = {"SIN2_THETA12": "SIN2_THETA12_SEED", "BETA_RAD": None,
               "DELTA_CKM_RAD": "DELTA_CKM_RAD"}
    ok_reg = True
    for o in obs:
        rid = reg_map.get(o["id"], o["id"])
        if rid is None:
            continue
        ok_reg &= abs(o["pred_val"] - float(reg[rid])) <= 1e-9 * abs(float(reg[rid]))
    check("scored predictions == frozen registry values (v84 lock, 13 obs)", ok_reg)

    # (A) exact census at the declared budgets
    cen = run_census(obs, G1_PQ_MAX, G2_PQ_MAX, DECADE)
    log10_joint = 0.0
    print("  --- per-observable hit densities (exact census) ---")
    for o in obs:
        h, d, p = cen[o["id"]]
        log10_joint += math.log10(p)
        print(f"    {o['id']:<16} window {o['window']:9.2e}  "
              f"p_i = {h:>9.0f}/{d:>11.0f} = {p:9.3e}   [{o['src']}]")
    print(f"  joint formula-fishing probability prod p_i = 1e{log10_joint:.2f}")

    check("fairness: every TFPT formula is inside its declared grammar class "
          "(13/13 reconstruct)", fairness_reconstruct(obs))
    check("every per-observable density p_i < 1 (the class never trivialises)",
          all(cen[o["id"]][2] < 1.0 for o in obs))
    check("joint null probability prod p_i < 1e-12 (look-elsewhere corrected)",
          log10_joint < -12.0)

    # (A') robustness: budget and window variations stay in the same regime
    var = []
    for pq1, pq2, dec in ((80, 8, 10.0), (200, 16, 10.0),
                          (130, 12, 5.0), (130, 12, 20.0)):
        c2 = run_census(obs, pq1, pq2, dec)
        var.append(sum(math.log10(c2[o["id"]][2]) for o in obs))
    spread = max(abs(v - log10_joint) for v in var)
    print(f"  robustness: log10 joint under budget/window variation: "
          f"{[f'{v:.1f}' for v in var]} (max shift {spread:.2f} dex)")
    check("census robust: coefficient bound 80..200 and window 5x..20x shift "
          "log10(prod p_i) by < 3 dex", spread < 3.0)

    # (B) Monte Carlo pseudo-theories (deterministic seeds)
    cA = run_mc(obs_by_id, MC_TRIALS, RNG_SEED, random_seed_constant=False)
    cB = run_mc(obs_by_id, MC_TRIALS, RNG_SEED + 1, random_seed_constant=True)
    maxA = max(k for k, n in enumerate(cA) if n)
    maxB = max(k for k, n in enumerate(cB) if n)
    meanA = sum(k * n for k, n in enumerate(cA)) / MC_TRIALS
    meanB = sum(k * n for k, n in enumerate(cB)) / MC_TRIALS
    print(f"  MC-A (phi0 fixed, random dials): mean K = {meanA:.3f}, "
          f"max K = {maxA}/13 over {MC_TRIALS} trials; dist {cA[:maxA+1]}")
    print(f"  MC-B (random seed phi*):         mean K = {meanB:.3f}, "
          f"max K = {maxB}/13 over {MC_TRIALS} trials; dist {cB[:maxB+1]}")
    check(f"MC-A: no random-dial theory reaches 10/13 hits in {MC_TRIALS} "
          f"trials (TFPT: 13/13); observed max {maxA}", maxA < 10)
    check(f"MC-B: no random-seed theory reaches 10/13 hits in {MC_TRIALS} "
          f"trials; observed max {maxB}", maxB < 10)
    check("MC means are far below the scorecard (mean K < 2 in both runs)",
          meanA < 2.0 and meanB < 2.0)

    # (C) negative controls -- the test has power
    tfpt_vals = [o["pred_val"] for o in
                 (obs_by_id[i] for i in MC_ORDER)]
    check("control: TFPT's own dial settings score 13/13 (plumbing)",
          hit_count(tfpt_vals, obs_by_id), 13, exact=True)

    def k_at(eps):
        best = 0
        for s in (+1, -1):
            f = PHI0 * (1.0 + s * eps)
            l = math.sqrt(f * (1.0 - f))
            vals = [obs_by_id[i]["pred"](f, l) for i in MC_ORDER]
            best = max(best, hit_count(vals, obs_by_id))
        return best
    k001, k01, k03 = k_at(0.01), k_at(0.10), k_at(0.30)
    print(f"  seed-perturbation collapse: K(+-1%)={k001}, K(+-10%)={k01}, "
          f"K(+-30%)={k03}  (vs 13/13 at phi0)")
    check("control: perturbing phi0 by 1% collapses the scorecard (K <= 9)",
          k001 <= 9)
    check("control: perturbing phi0 by 10% collapses the scorecard (K <= 6)",
          k01 <= 6)

    rng = random.Random(RNG_SEED + 2)
    tot = mx = 0
    for _ in range(1000):
        perm = list(range(13))
        rng.shuffle(perm)
        k = sum(1 for i, j in enumerate(perm)
                if abs(tfpt_vals[i] - obs_by_id[MC_ORDER[j]]["m"])
                <= obs_by_id[MC_ORDER[j]]["window"] * obs_by_id[MC_ORDER[j]]["m"])
        tot += k
        mx = max(mx, k)
    print(f"  shuffled-data control: mean K = {tot/1000:.2f}, max K = {mx} "
          f"over 1000 permutations")
    check("control: shuffling data<->formula assignment collapses the "
          "scorecard (mean K <= 2, max K <= 6)", tot / 1000 <= 2.0 and mx <= 6)

    # (D) the alpha cubic class
    n_alpha, hits_alpha, tfpt_ainv = alpha_class_census()
    p_alpha = max(hits_alpha, 1) / n_alpha
    print(f"  alpha class: {n_alpha} F_U1 variants root-solved; "
          f"{hits_alpha} inside the 4e-8 CODATA window; p_alpha <= {p_alpha:.1e}")
    check("alpha fairness: the TFPT F_U1 instance is in the class and "
          "reproduces the frozen root", abs(tfpt_ainv - 137.0359992168407)
          < 1e-6)
    check("alpha look-elsewhere: at most 3 of the F_U1 variants hit the "
          f"4e-8 CODATA window (observed {hits_alpha})",
          1 <= hits_alpha <= 3)
    check("alpha null probability p_alpha < 1e-4", p_alpha < 1e-4)

    # headline
    log10_total = log10_joint + math.log10(p_alpha)
    print(f"  HEADLINE: P(random complexity-matched theory reproduces the "
          f"scorecard incl. alpha) <= 1e{log10_total:.1f}")
    print(f"            = {-log10_total * math.log2(10):.0f} bits of "
          f"look-elsewhere-corrected surprise (conditional on the declared "
          f"grammar; never 'certainty')")
    check("headline: joint null probability (13 obs + alpha) < 1e-16",
          log10_total < -16.0)

    return summary("v100 numerology null test")


def fairness_reconstruct(obs):
    """Re-evaluate every TFPT formula from its declared grammar encoding."""
    ok = True
    for o in obs:
        if o["arity"] == 1:
            p, q, a, c, d, b, u, w = o["g1"]
            ok &= (1 <= p <= G1_PQ_MAX and 1 <= q <= G1_PQ_MAX
                   and a in G1_A and c in G1_C and d in G1_D and b in G1_B
                   and abs(u) <= G1_EXP_NUM and 1 <= w <= G1_EXP_DEN)
            v = (p / q) * PHI0 ** a * LAM ** c * (1 + LAM) ** d * PI ** b \
                * math.exp(u / w)
        else:
            (p1, q1, a1, c1, b1), sgn, (p2, q2, a2, c2, b2) = o["g2"]
            ok &= all(1 <= x <= G2_PQ_MAX for x in (p1, q1, p2, q2))
            ok &= (a1 in G2_A and c1 in G2_C and b1 in G2_B
                   and a2 in G2_A and c2 in G2_C and b2 in G2_B)
            v = ((p1 / q1) * PHI0 ** a1 * LAM ** c1 * PI ** b1
                 + sgn * (p2 / q2) * PHI0 ** a2 * LAM ** c2 * PI ** b2)
        ok &= abs(v - o["pred_val"]) <= 1e-12 * max(abs(o["pred_val"]), 1e-30)
    return ok


if __name__ == "__main__":
    raise SystemExit(1 if run() else 0)

8eeee76c7fc3789f8ce7777ff62835d2cf4aa37dc2e28cb5850480faf72f94f5

The public registry includes 116 document views, 6,693 section records and 985 physical LaTeX sources. Active publication documents, research material and archives retain their recorded roles. This is source navigation, not a new semantic review or proof run.