# build_math_page.py - writes build/math.html ("The math, step by step").
#
# NOTHING on that page is typed.  Every constant and every measured reference value is read from math_inputs.json (each with its
# source), every derived value is computed here (logic ported from downloads/code/*.py), and every PASS / FAIL / PARTIAL badge is
# decided by code from a tolerance that is stated on the page.  After writing, a provenance check scans the page and refuses
# (non-zero exit) if any number is not one that this script emitted or part of an equation's own formula.
#
#   python -I build_math_page.py <build_dir>              # build + provenance check
#   python -I build_math_page.py <build_dir> --plant      # also type a number into a copy and prove the check catches it
import sys
sys.dont_write_bytecode = True
import html, importlib.util, json, math, pathlib, re, zipfile
import xml.etree.ElementTree as ET
import numpy as np
from scipy.integrate import quad, solve_ivp
from scipy.optimize import brentq

HERE = pathlib.Path(__file__).resolve().parent
spec = importlib.util.spec_from_file_location("site_common", HERE / "site_common.py")
sc = importlib.util.module_from_spec(spec); spec.loader.exec_module(sc)
BUILD = pathlib.Path(sys.argv[1])
INP = json.loads((HERE / "math_inputs.json").read_text(encoding="utf-8"))
def _walk(o):
    if isinstance(o, str):
        yield o
    elif isinstance(o, dict):
        for v in o.values(): yield from _walk(v)
    elif isinstance(o, list):
        for v in o: yield from _walk(v)
LOADED_TEXT = " ".join(_walk(INP))
def V(k): return INP[k]["value"]
def CL(k): return INP["record_claims"]["value"][k]
STDM = V("std_ns_mass")
def SRC(k): return INP[k]["source"]

# ------------------------------------------------------------------------------------------------ static check of THIS file
import ast
def static_check():
    """No string literal that reaches the page may contain a typed number, and no N()/I()/fr() call may take a numeric literal,
    except: equation strings (the 5th argument of check()), the ANCH anchor table, and the module docstring."""
    src = pathlib.Path(__file__).read_text(encoding="utf-8"); tree = ast.parse(src); bad = []
    exempt = set()
    for node in ast.walk(tree):
        if isinstance(node, ast.Call) and getattr(node.func, "id", "") == "check" and len(node.args) > 4 and isinstance(node.args[4], ast.List):
            for d in ast.walk(node.args[4]): exempt.add(id(d))
        if isinstance(node, ast.Assign) and any(getattr(tg, "id", "") in ("ANCH", "NUM", "SUP") for tg in node.targets):
            for d in ast.walk(node): exempt.add(id(d))
        if isinstance(node, ast.Call) and getattr(node.func, "id", "") == "eq":
            for d in ast.walk(node): exempt.add(id(d))
    for node in ast.walk(tree):
        if id(node) in exempt: continue
        if isinstance(node, ast.Constant) and isinstance(node.value, str) and typed_numbers(node.value):
            if node.value.startswith("# build_math_page") : continue
            bad.append((node.lineno, node.value[:60]))
        if isinstance(node, ast.Call) and getattr(node.func, "id", "") in ("N", "I") and node.args and isinstance(node.args[0], ast.Constant) and isinstance(node.args[0].value, (int, float)):
            bad.append((node.lineno, "numeric literal passed to " + node.func.id))
    return bad
def typed_numbers(s):
    """decimals and numbers of three or more digits typed into a string, unless the loaded data file itself contains them."""
    loaded = set(NUM.findall(LOADED_TEXT))
    return [tk for tk in NUM_S.findall(s) if ("." in tk or len(tk) >= 3) and tk not in loaded]
NUM_S = re.compile(r"(?<![A-Za-z_\d.\-])\d+(?:\.\d+)?")

# ------------------------------------------------------------------------------------------------ provenance machinery
NUM = re.compile(r"(?<![A-Za-z_\d.])\d+(?:\.\d+)?")
EMITTED, FORMULA = set(NUM.findall(LOADED_TEXT)), set()
def emit(s):
    EMITTED.update(NUM.findall(s)); return s
def eq(s):
    FORMULA.update(NUM.findall(s)); return s
SUP = str.maketrans("0123456789-", "⁰¹²³⁴⁵⁶⁷⁸⁹⁻")
def N(x, sig=4):
    """format a computed or loaded number; every digit group goes into the emitted set."""
    x = float(x)
    if x == 0: return emit("0")
    sign = "−" if x < 0 else ""; a = abs(x)
    if 1e-3 <= a < 1e6:
        dec = max(0, sig - 1 - int(math.floor(math.log10(a))))
        s = f"{a:.{dec}f}"
        if "." in s: s = s.rstrip("0").rstrip(".")
    else:
        e = int(math.floor(math.log10(a))); m = a / 10 ** e
        if round(m, sig - 1) >= 10: m /= 10; e += 1
        ms = f"{m:.{sig - 1}f}"
        if "." in ms: ms = ms.rstrip("0").rstrip(".")
        s = ms + "×10" + str(e).translate(SUP)
    return emit(sign + s)
def I(x): return emit(str(int(x)))
def fr(a, b): return emit(f"{int(a)}/{int(b)}")

if "--static" in sys.argv:
    b = static_check()
    for ln, s in b: print(ln, ascii(s))
    print("static offenders:", len(b)); sys.exit(1 if b else 0)
STATUSES = ["PASS", "FAIL", "PARTIAL", "NOT TESTABLE YET", "INPUT", "INHERITED", "NOT RE-RUN"]
CLS = {"PASS": "pass", "FAIL": "fail", "PARTIAL": "partial", "NOT TESTABLE YET": "nty", "INPUT": "input", "INHERITED": "inherited", "NOT RE-RUN": "norerun"}
PF = V("partial_factor")
def j_rel(x, ref, tol):
    d = abs(x - ref) / abs(ref)
    return "PASS" if d <= tol else ("PARTIAL" if d <= PF * tol else "FAIL")
def j_sigma(x, ref, sig, n=1):
    d = abs(x - ref) / sig
    return "PASS" if d <= n else ("PARTIAL" if d <= PF * n else "FAIL")
def j_range(x, lo, hi): return "PASS" if lo <= x <= hi else "FAIL"
MISMATCH = []
def rec_line(label, computed, claimed, sig=4):
    """the record's own stated number beside this run's number; a >1% difference is flagged and counted."""
    d = abs(computed - claimed) / abs(claimed) if claimed else abs(computed)
    ok = d <= V("tol_percent_1") or round(computed, 3) == claimed
    if not ok: MISMATCH.append(label)
    return (f"The record states {N(claimed, sig)}; this run gives {N(computed, sig)}. " + ("They agree." if ok else "They DIFFER."))

# ------------------------------------------------------------------------------------------------ record anchors
R = "record.html#"
ANCH = {
 "A": "item-a1", "K1": "item-k1", "K2": "item-k2", "K3": "item-k3", "K4": "item-k4", "K5": "item-k5",
 "K6": "k6.-each-particles-own-scale-view-matthew-13-nov", "K7": "k7.-the-particle-mass-formula-fixed-for-the-leptons",
 "K8": "k8.-geminis-other-proton-results-checked-line-by-line", "K9": "k9.-black-hole-ringdown-and-the-central-particle",
 "K10": "k10.-antiparticle-matrices-and-decay-rates-matthews-saved-requests-workbook-updated",
 "M1": "m1.-the-containment-potential-vp-pinned-with-no-knobs",
 "M2": "m2.-the-mass-couplings-y-written-as-relations-between-generations",
 "M3": "m3.-the-layer-structure-g_ab-matthews-s-is-the-only-fair-distance",
 "M5": "m5.-space-as-crossing-waves-why-the-geometry-keeps-transitioning-matthew-original-idea-newly-written-down-14-nov",
 "M6": "m6.-quarks-a-hint-not-yet-a-result",
 "M7": "m7.-the-cause-of-the-waves-energy-spreading-through-uneven-resistance-matthew-14-nov",
 "M8": "m8.-the-cosmos-test-t4-first-pass", "M9": "m9.-the-factor-2-what-kind-of-resistance-gravity-must-be",
 "M10": "m10.-gravitys-other-two-properties-reaching-everywhere-and-acting-on-everything-equally",
 "M11": "m11.-the-expansion-equation-the-full-bag-and-the-full-glass-matthew-original-theory",
 "M12": "m12.-entanglement-and-uncertainty-inside-the-bag",
 "M14": "m14.-gravity-as-a-byproduct-of-entanglement-matthew-14-nov-gravity-is-accidental-as-a-result-of-entanglement",
 "M15": "m15.-quantum-physics-is-entanglement-the-universe-as-one-connected-web-matthew-from-his-original-meditations",
 "M16": "m16.-entangled-particles-are-touching-the-area-law-charges-from-consistency",
 "M17": "m17.-final-scorecard-14-nov-2026", "M18": "m18.-the-last-two-open-problems-honest-best-attempts",
 "P1": "p1.-the-master-equation-written-out-in-full",
 "Q3": "q3.-b-set-by-the-proton-with-nothing-adjusted-script-math_fixesproton_sets_b_2026-11-14.py",
 "Q4": "q4.-neutron-stars-with-nothing-adjusted-script-math_fixesneutron_stars_no_knobs_2026-11-14.py",
 "Q5": "q5.-quark-masses-one-idea-tested-and-it-fails-script-math_fixesquark_link_pattern_2026-11-14.py",
 "P2": "p2.-quark-masses-at-one-common-energy-script-math_fixesquark_common_scale_2026-11-14.py",
 "P4": "p4.-gravity-tightened-script-math_fixesgravity_tighten_2026-11-14.py",
 "P5": "p5.-the-29-twist-and-the-12-link-script-math_fixestwist_and_link_2026-11-14.py",
}

CHECKS = []
SECTIONS = {}
def check(cid, sec, title, rec, eqs, subst, computed, measured, tol, status, note="", pre=None):
    if pre: emit(pre)
    assert status in STATUSES, status
    CHECKS.append(dict(id=cid, sec=sec, title=title, rec=rec, eqs=[eq(e) for e in eqs], subst=subst, computed=computed,
                       measured=measured, tol=tol, status=status, note=note, pre=pre))

# ------------------------------------------------------------------------------------------------ computations
G, C = V("G_SI"), V("c_SI"); MSUN, RSUN = V("Msun_kg"), V("Rsun_m"); ASEC = V("arcsec_per_rad")
me, mmu, mtau = V("m_e"), V("m_mu"), V("m_tau")
Q23 = V("koide_target_num") / V("koide_target_den"); D29 = V("brannen_delta_num") / V("brannen_delta_den")
def Qk(a, b, c): return (a + b + c) / (math.sqrt(a) + math.sqrt(b) + math.sqrt(c)) ** 2

def bending():
    b = RSUN; z = np.linspace(-2e12, 2e12, 2_000_001); r = np.sqrt(b * b + z * z)
    dndb = -2 * G * MSUN / (C * C) * b / r ** 3
    return -np.trapezoid(dndb, z) * ASEC

def age_gyr():
    Om, Or = V("Omega_m"), V("Omega_r"); OL = 1 - Om - Or
    H = lambda a: V("H0") * math.sqrt(Om / a ** 3 + Or / a ** 4 + OL)
    return quad(lambda a: 1 / (a * H(a)), 1e-10, 1)[0] * V("Mpc_km") / V("Gyr_s"), OL, H

# --- TOV (units: p, e in MeV/fm^3, r and m in km) ---
CONV, MSUN_KM = V("tov_conv"), V("Msun_geo_km")
def tov(pf, einv, ec, psurf=1e-6):
    pc = pf(ec)
    def rhs(r, y):
        p, m = y
        if p <= psurf: return [0, 0]
        e = einv(p) * CONV; pg = p * CONV
        return [-(e + pg) * (m + 4 * np.pi * r ** 3 * pg) / (r * (r - 2 * m)) / CONV, 4 * np.pi * r * r * e]
    ev = lambda r, y: y[0] - psurf; ev.terminal = True
    r0 = 1e-4
    s = solve_ivp(rhs, [r0, 60], [pc, 4 / 3 * np.pi * r0 ** 3 * ec * CONV], events=ev, rtol=1e-8, atol=1e-12, max_step=0.02)
    return s.t[-1], s.y[1][-1] / MSUN_KM

def bag_eos(B, av):
    pf = lambda e: max((e - 4 * B) / 3 + av * (e - 4 * B) ** 2 / (1 + 1.45 * av * (e - 4 * B)), 0.0)
    def einv(p):
        if p <= 0: return 4 * B
        if av == 0: return 3 * p + 4 * B
        A = 1.45 * av / 3 + av; Bc = 1 / 3 - 1.45 * av * p
        return 4 * B + (-Bc + math.sqrt(Bc * Bc + 4 * A * p)) / (2 * A)
    return pf, einv

def bag_sequence(B, av):
    pf, einv = bag_eos(B, av)
    rows = np.array([(ec,) + tov(pf, einv, ec, psurf=1e-7) for ec in np.linspace(4 * B * 1.01, V("tov_scan_max"), 60)])
    i = rows[:, 2].argmax(); r = rows[:i + 1]
    return float(rows[i, 2]), float(np.interp(STDM, r[:, 2], r[:, 1])), float(rows[i, 1])

def gemini_sequence():
    B, a = V("bag_B_gemini"), V("bag_wall_gemini"); sc_, pw = V("gemini_p_scale"), V("gemini_p_power")
    def pf(e):
        if e <= 4 * B: return max(1e-4, sc_ * e * (e / (4 * B)) ** pw)
        d = e - 4 * B; return d / 3 + a * d * d / (1 + 1.45 * a * d)
    egrid = np.logspace(-1, np.log10(5000), 40000)
    pg = np.maximum.accumulate(np.array([pf(e) for e in egrid]) + egrid * 1e-12)
    einv = lambda p: np.interp(p, pg, egrid)
    rows = np.array([(ec,) + tov(pf, einv, ec) for ec in sorted(V("gemini_ns_ec_list_low") + V("gemini_ns_ec_list"))])
    i = rows[:, 2].argmax(); r = rows[:i + 1]
    return float(rows[i, 2]), float(np.interp(STDM, r[:, 2], r[:, 1])), pf(150.0)

def sly_tables():
    c2 = V("c2_cgs"); crust = V("sly_crust"); G1, G2, G3 = V("sly_G"); lr1, lr2 = V("sly_log_rho12")
    rho1, rho2 = 10 ** lr1, 10 ** lr2; p1 = 10 ** V("sly_log_p1") / c2; K1 = p1 / rho1 ** G1
    K4c, G4c = crust[3][:2]; rho0 = (K4c / K1) ** (1 / (G1 - G4c))
    pieces = []; lo = 0
    for K, Gm, hi in crust[:3]: pieces.append([K, Gm, lo, hi]); lo = hi
    pieces.append([K4c, G4c, lo, rho0]); pieces.append([K1, G1, rho0, rho1])
    K2 = K1 * rho1 ** (G1 - G2); pieces.append([K2, G2, rho1, rho2])
    K3 = K2 * rho2 ** (G2 - G3); pieces.append([K3, G3, rho2, 1e17])
    prev = None
    for pc in pieces:
        K, Gm, lo, hi = pc
        if prev is None: pc.append(0.0)
        else:
            Kp, Gp, _, _, ap = prev; e_lo = (1 + ap) * lo + Kp / (Gp - 1) * lo ** Gp
            pc.append(e_lo / lo - 1 - K / (Gm - 1) * lo ** (Gm - 1))
        prev = pc
    conv = V("mevfm3_gcm3"); rr = np.logspace(3, 16.5, 60000); P = []; E = []
    for r in rr:
        for K, Gm, lo, hi, ai in pieces:
            if lo <= r < hi: P.append(K * r ** Gm); E.append((1 + ai) * r + K / (Gm - 1) * r ** Gm); break
    return np.array(E) / conv, np.array(P) / conv   # note: arrays are built so that (E, P) rise together

def k2_of(Cc, y):
    t1 = 8 * Cc ** 5 / 5 * (1 - 2 * Cc) ** 2 * (2 + 2 * Cc * (y - 1) - y)
    t2 = (2 * Cc * (6 - 3 * y + 3 * Cc * (5 * y - 8)) + 4 * Cc ** 3 * (13 - 11 * y + Cc * (3 * y - 2) + 2 * Cc * Cc * (1 + y))
          + 3 * (1 - 2 * Cc) ** 2 * (2 - y + 2 * Cc * (y - 1)) * np.log(1 - 2 * Cc))
    return t1 / t2

def tidal_star(pf, einv, dpde, ec, es_surf, psurf=1e-9):
    pc = pf(ec)
    def rhs(r, Y):
        p, m, y = Y
        e_ = einv(max(p, psurf)); e = e_ * CONV; P = p * CONV
        f = 1 - 2 * m / r; cs2 = dpde(e_)
        dp = -(e + P) * (m + 4 * np.pi * r ** 3 * P) / (r * r * f)
        Q = 4 * np.pi / f * (5 * e + 9 * P + (e + P) / cs2) - 6 / (r * r * f) - (2 * (m + 4 * np.pi * r ** 3 * P) / (r * r * f)) ** 2
        F = (1 - 4 * np.pi * r * r * (e - P)) / f
        return [dp / CONV, 4 * np.pi * r * r * e, -(y * y + y * F + r * r * Q) / r]
    ev = lambda r, Y: Y[0] - psurf; ev.terminal = True
    r0 = 1e-4
    s = solve_ivp(rhs, [r0, 60], [pc, 4 / 3 * np.pi * r0 ** 3 * ec * CONV, 2.0], events=ev, rtol=1e-9, atol=1e-14, max_step=0.01)
    R_ = s.t[-1]; M = s.y[1][-1]; y = s.y[2][-1]
    y -= 4 * np.pi * R_ ** 3 * es_surf * CONV / M
    Cc = M / R_; k2 = k2_of(Cc, y)
    return M / MSUN_KM, R_, k2, 2 / 3 * k2 / Cc ** 5

def bag_tidal(B, av, n=45):
    pf, einv = bag_eos(B, av); h = 1e-3
    dpde = lambda e: (pf(e + h) - pf(max(e - h, 4 * B))) / (e + h - max(e - h, 4 * B))
    rows = np.array([tidal_star(pf, einv, dpde, ec, 4 * B) for ec in np.linspace(4 * B * 1.05, 2500, n)])
    i = rows[:, 0].argmax(); r = rows[:i + 1]
    return float(rows[i, 0]), float(np.interp(STDM, r[:, 0], r[:, 3])), float(np.interp(STDM, r[:, 0], r[:, 2]))

# --- minimal xlsx reader (no pandas in the isolated interpreter) ---
def read_xlsx(path, sheet):
    ns = {"m": "http://schemas.openxmlformats.org/spreadsheetml/2006/main"}
    with zipfile.ZipFile(path) as z:
        ss = []
        if "xl/sharedStrings.xml" in z.namelist():
            for si in ET.fromstring(z.read("xl/sharedStrings.xml")).findall("m:si", ns):
                ss.append("".join(t.text or "" for t in si.iter("{%s}t" % ns["m"])))
        wb = ET.fromstring(z.read("xl/workbook.xml"))
        rels = {r.get("Id"): r.get("Target") for r in ET.fromstring(z.read("xl/_rels/workbook.xml.rels"))}
        rid = None
        for s in wb.find("m:sheets", ns):
            if s.get("name") == sheet: rid = s.get("{http://schemas.openxmlformats.org/officeDocument/2006/relationships}id")
        target = rels[rid].lstrip("/"); target = target if target.startswith("xl/") else "xl/" + target
        rows = []
        for row in ET.fromstring(z.read(target)).find("m:sheetData", ns):
            vals = {}
            for c in row:
                col = re.match(r"[A-Z]+", c.get("r")).group(0); t = c.get("t"); v = c.find("m:v", ns)
                if t == "inlineStr": val = "".join(x.text or "" for x in c.iter("{%s}t" % ns["m"]))
                elif v is None: val = None
                elif t == "s": val = ss[int(v.text)]
                elif t in ("str", "b", "e"): val = v.text
                else: val = float(v.text)
                vals[col] = val
            rows.append(vals)
    head = rows[0]; cols = sorted(head, key=lambda k: (len(k), k))
    return [{head[c]: r.get(c) for c in cols} for r in rows[1:]]

XLSX = BUILD / "downloads" / "TRPG_Every_Particle_to_Every_Particle.xlsx"
# ---- rev3: running the published scripts as published (python -I, separate working directory) and the no-knobs neutron-star port
import subprocess, shutil
from scipy.interpolate import interp1d
def arg_after(flag): return sys.argv[sys.argv.index(flag) + 1] if flag in sys.argv else None
OUTDIR = pathlib.Path(arg_after("--outputs")) if arg_after("--outputs") else None
PUBLISHED = ("proton_sets_B", "neutron_stars_no_knobs", "quark_link_pattern", "area_law_and_charges", "expansion_and_entanglement",
             "master_action_M1_containment_potential", "neutron_star_tidal_fix", "build_antiparticle_and_decays")
def clean_paths(s):
    """no local file-system paths on the public page: keep only the last component of any Windows path"""
    return re.sub(r"[A-Za-z]:[\\/][^\s\"']+", lambda m: "..." + "\\" + m.group(0).replace("/", "\\").split("\\")[-1], s)
def published(prefix):
    f_, rc_, o_, e_ = published_raw(prefix); return f_, rc_, clean_paths(o_), clean_paths(e_)
def published_raw(prefix):
    f = sorted((BUILD / "downloads" / "code").glob(prefix + "*.py"))[0]
    if OUTDIR is not None and (OUTDIR / (f.stem + ".txt")).exists():
        rc = None
        for ln in (OUTDIR / "_done.txt").read_text(encoding="utf-8").splitlines():
            if ln.startswith(f.stem + " rc="): rc = int(ln.split("rc=")[1].split()[0])
        return f.name, rc, (OUTDIR / (f.stem + ".txt")).read_text(encoding="utf-8", errors="replace"), (OUTDIR / (f.stem + ".err")).read_text(encoding="utf-8", errors="replace")
    run = HERE / V("script_rundir"); run.mkdir(exist_ok=True)
    shutil.copyfile(BUILD / "downloads" / "TRPG_Every_Particle_to_Every_Particle.xlsx", run / "TRPG_Every_Particle_to_Every_Particle.xlsx")
    pr = subprocess.run([sys.executable, "-I", str(f)], cwd=str(run), capture_output=True, text=True, timeout=7200)
    return f.name, pr.returncode, pr.stdout, pr.stderr
PUB = {p: published(p) for p in PUBLISHED}
for _p in PUB.values(): emit(_p[2]); emit(_p[3][-1500:])
def nums(line): return [float(x) for x in re.findall(r"-?\d+\.\d+|-?\d+", line)]

_nkf = sorted((BUILD / "downloads" / "code").glob("neutron_stars_no_knobs*.py"))[0]
_nkspec = importlib.util.spec_from_file_location("nk_published", _nkf); NKMOD = importlib.util.module_from_spec(_nkspec); _nkspec.loader.exec_module(NKMOD)
def nk_sequence(B, X, gap):
    """Matt's own sequence() from the published neutron_stars_no_knobs script, called directly (its main block is not run)."""
    M_, L_, R_, a4_, es_ = NKMOD.sequence(B, X, gap)
    return float(M_), float(L_), float(R_)


# ================================================================================================ SECTION A: home-page table
SEC_A = "The home-page results table and the neutrino prediction"
SECTIONS["A"] = SEC_A
bend = bending()
check("A1", "A", "Bending of starlight by the Sun", ("results.html", "Results page"),
      ["δ = ∫ ∇⊥n dz,  with n = 1 + 2GM/(rc²)"],
      f"M = {N(MSUN)} kg, R = {N(RSUN)} m, G = {N(G)}, c = {N(C)} → n − 1 = {N(2*G*MSUN/(RSUN*C*C))} at the Sun's edge",
      [f"Deflection = {N(bend)} arcsec (numerical integral along the ray)", rec_line("A1 bending", bend, CL("A1"))],
      f"{N(V('bending_measured'))} arcsec — {SRC('bending_measured')}", f"within {N(V('tol_percent_1')*100,2)}% of the measured value",
      j_rel(bend, V("bending_measured"), V("tol_percent_1")))
age, OLs, Hf = age_gyr()
check("A2", "A", "Age of the universe", (R + ANCH["M11"], "record M11"),
      ["t = ∫ da / (a H(a)),  H = H₀ √(Ω_m a⁻³ + Ω_r a⁻⁴ + Ω_Λ)"],
      f"H₀ = {N(V('H0'))} km/s/Mpc, Ω_m = {N(V('Omega_m'))}, Ω_r = {N(V('Omega_r'))}, Ω_Λ = 1 − Ω_m − Ω_r = {N(OLs,5)}",
      [f"Age = {N(age)} Gyr", rec_line("A2 age", age, CL("A2"))],
      f"{N(V('age_measured'))} ± {N(V('age_sigma'))} Gyr — {SRC('age_measured')}", "within 1 standard deviation of the measurement",
      j_sigma(age, V("age_measured"), V("age_sigma")))
bB = V("bag_B_gemini")
mm, r14, _ = bag_sequence(bB, 0.0)
# ---- rev3 A3: neutron stars with nothing adjusted (record Q4, home page)
Bp = V("bag_B_proton"); gaps = V("pairing_gap_range"); Xc = V("no_knobs_scale")
nk = {(Bv, gp): nk_sequence(Bv, Xc, gp) for Bv in Bp for gp in gaps}
nkM = [v[0] for v in nk.values()]; nkL = [v[1] for v in nk.values()]; nkR = [v[2] for v in nk.values()]
nk_sub = (f"B = {', '.join(N(b) for b in Bp)} MeV/fm\u00b3 (set by the proton, stage B9); pairing gap \u0394 = {N(gaps[0])} and {N(gaps[1])} MeV; scale factor X = {N(Xc)}; "
          f"strange quark mass {N(V('no_knobs_ms'))} MeV; \u03b1_s from {N(V('no_knobs_alpha_s_MZ'))} at M_Z by two-loop running")
nk_eq = ["p(\u03bc) = (3/4\u03c0\u00b2)[a\u2084\u03bc\u2074 \u2212 (m_s\u00b2 \u2212 4\u0394\u00b2)\u03bc\u00b2] \u2212 B,   a\u2084 = 1 \u2212 2\u03b1_s/\u03c0;   then the TOV equations"]
mt, ms_ = V("ns_mass_test"), V("ns_mass_sigma")
check("A3a", "A", "Heaviest neutron star (B from the proton, quark pairing)", (R + ANCH["Q4"], "record Q4"), nk_eq, nk_sub,
      [f"Maximum mass from {N(min(nkM),3)} to {N(max(nkM),3)} solar masses"] + [f"B = {N(k_[0])}, Δ = {N(k_[1])} MeV: M_max = {N(v_[0],3)}, Λ = {N(v_[1],3)}, R = {N(v_[2],3)} km" for k_, v_ in nk.items()] + [rec_line("A3a low", min(nkM), CL("Q4mlo"), 3), rec_line("A3a high", max(nkM), CL("Q4mhi"), 3)],
      f"at least {N(mt)} solar masses (the measured star is {N(V('ns_mass_measured'))} \u00b1 {N(ms_)}) \u2014 {SRC('ns_mass_test')}", "every case in the proton band at or above the test mass",
      "PASS" if min(nkM) >= mt else ("PARTIAL" if min(nkM) >= mt - ms_ else "FAIL"))
rlo, rhi = V("ns_radius_lo"), V("ns_radius_hi")
inside = [rlo <= r_ <= rhi for r_ in nkR]
check("A3b", "A", "Radius of a standard-mass neutron star (same cases)", (R + ANCH["Q4"], "record Q4"), ["same star sequences; R at M = 1.4 solar masses"], nk_sub,
      [f"R({N(STDM)}) from {N(min(nkR),3)} to {N(max(nkR),3)} km", rec_line("A3b low", min(nkR), CL("Q4Rlo"), 3), rec_line("A3b high", max(nkR), CL("Q4Rhi"), 3)],
      f"{N(rlo)} to {N(rhi)} km \u2014 {SRC('ns_radius_lo')}", "every case inside the measured range",
      "PASS" if all(inside) else ("PARTIAL" if any(inside) else "FAIL"),
      "Part of the predicted range sits below the lower edge of the quoted measurement." if not all(inside) else "")
check("A3c", "A", "Tidal deformability \u039b of a standard-mass neutron star (same cases)", (R + ANCH["Q4"], "record Q4"), ["\u039b = (2/3) k\u2082 / C\u2075, from the same sequences (log-interpolated in \u039b)"], nk_sub,
      [f"\u039b({N(STDM)}) from {N(min(nkL),3)} to {N(max(nkL),3)}", rec_line("A3c low", min(nkL), CL("Q4Llo"), 3), rec_line("A3c high", max(nkL), CL("Q4Lhi"), 3)],
      f"\u039b at or below {N(V('lambda_limit'))} \u2014 {SRC('lambda_limit')}", "every case at or below the limit", "PASS" if max(nkL) <= V("lambda_limit") else "FAIL")
rows_nk = [nums(ln) for ln in PUB["neutron_stars_no_knobs"][2].splitlines() if re.match(r"^\s+\d+\.\d\s+\d\.\d\s+\d\.\d{3}\s", ln)]
strong = [r_[3] for r_ in rows_nk[:9]] if len(rows_nk) >= 9 else []
check("A3d", "A", "The strong force alone, with no quark pairing", (R + ANCH["Q4"], "record Q4"), ["same equation of state with \u0394 = 0"],
      "from the published script's own printed table (nine cases: three values of B, three scale factors)",
      ([f"Maximum mass from {N(min(strong),3)} to {N(max(strong),3)} solar masses", rec_line("A3d low", min(strong), CL("Q4strong_lo"), 3), rec_line("A3d high", max(strong), CL("Q4strong_hi"), 3)] if strong else ["The published script's table could not be read."]),
      f"at least {N(mt)} solar masses \u2014 {SRC('ns_mass_test')}", "heaviest star at or above the test mass", ("PASS" if strong and min(strong) >= mt else "FAIL"),
      "Fails the two-solar-mass test, as the record says: this is why the pairing gap is needed.")
check("A3e", "A", "The predicted quark pairing gap", (R + ANCH["Q4"], "record Q4"), ["\u0394 such that both neutron-star tests pass"],
      f"gap range {N(gaps[0])} to {N(gaps[1])} MeV; theory estimates {N(V('pairing_theory_range')[0])} to {N(V('pairing_theory_range')[1])} MeV",
      ["No measurement of the pairing gap exists."], "not measured", "a measured pairing gap", "NOT TESTABLE YET",
      "A prediction that can fail: it has not been measured, as the record says.")
# leptons
Qlep = Qk(me, mmu, mtau)
a_, b_ = math.sqrt(me), math.sqrt(mmu); s_ = a_ + b_
A_, B_, C_ = 1 / 3, -4 / 3 * s_, (a_ * a_ + b_ * b_) - 2 / 3 * s_ * s_
x_ = (-B_ + math.sqrt(B_ * B_ - 4 * A_ * C_)) / (2 * A_); tau_pred = x_ * x_
Mt = math.sqrt(mtau) / (1 + math.sqrt(2) * math.cos(D29))
mu_pred = (Mt * (1 + math.sqrt(2) * math.cos(D29 + 4 * math.pi / 3))) ** 2
e_pred = (Mt * (1 + math.sqrt(2) * math.cos(D29 + 2 * math.pi / 3))) ** 2
tl = V("tol_lepton_mass")
check("A4a", "A", "Tau mass from the electron and muon (Koide)", (R + ANCH["K7"], "record K7"),
      ["(m_e + m_μ + m_τ) / (√m_e + √m_μ + √m_τ)² = 2/3, solved for m_τ"],
      f"m_e = {N(me,8)} MeV, m_μ = {N(mmu,10)} MeV; Q on the measured trio = {N(Qlep,7)}",
      [f"Predicted m_τ = {N(tau_pred,6)} MeV", rec_line("A4a tau", tau_pred, CL("A4a"), 6)],
      f"{N(mtau,6)} MeV — {SRC('m_tau')}", f"within {N(tl*100,2)}%", j_rel(tau_pred, mtau, tl))
check("A4b", "A", "Electron and muon masses from the tau alone", (R + ANCH["K7"], "record K7"),
      ["√m_k = M(1 + √2 cos(δ + 2πk/3)),  δ = 2/9,  M fixed by m_τ"],
      f"δ = {fr(V('brannen_delta_num'), V('brannen_delta_den'))} = {N(D29,6)}; M from m_τ = {N(mtau,6)} MeV",
      [f"Predicted m_e = {N(e_pred,6)} MeV ({N((e_pred/me-1)*100,2)}% off); predicted m_μ = {N(mu_pred,7)} MeV ({N((mu_pred/mmu-1)*100,2)}% off)",
       "Three masses from one input scale (plus the fixed twist)."],
      f"m_e = {N(me,8)} MeV, m_μ = {N(mmu,10)} MeV — {SRC('m_e')}", f"each within {N(tl*100,2)}%",
      "PASS" if j_rel(e_pred, me, tl) == "PASS" and j_rel(mu_pred, mmu, tl) == "PASS" else ("PARTIAL" if max(abs(e_pred/me-1), abs(mu_pred/mmu-1)) <= PF*tl else "FAIL"))
cab = math.sqrt(V("q_d") / V("q_s"))
check("A5", "A", "Cabibbo angle from quark masses (label: already known, Gatto, Sartori and Tonin 1968)", (R + ANCH["K4"], "record K4"),
      ["θ_C = e^S = √(m_d / m_s),  S = ½ ln(m_d / m_s)"],
      f"m_d = {N(V('q_d'))} MeV, m_s = {N(V('q_s'))} MeV ({SRC('q_d')})",
      [f"√(m_d/m_s) = {N(cab,4)}", rec_line("A5 cabibbo", cab, CL("A5"))],
      f"{N(V('ckm_Vus'))} — {SRC('ckm_Vus')}", f"within {N(V('tol_percent_1')*100,2)}% (the quark-mass uncertainty)", j_rel(cab, V("ckm_Vus"), V("tol_percent_1")),
      "The relation is Gatto, Sartori and Tonin (1968); it is known physics, not new here.")
# tau lifetime
GF, hb = V("GF_MeV2"), V("hbar_MeVs")
tau_mu_c = hb / (GF ** 2 * mmu ** 5 / (192 * math.pi ** 3)); S_tm = 0.5 * math.log(mtau / mmu)
tau_tau_c = tau_mu_c * V("BR_tau_e") * math.exp(-10 * S_tm) * 1e15
check("A6", "A", "Tau lifetime from the muon (label: already known, lepton universality)", (R + ANCH["K10"], "record K10"),
      ["τ_μ = 192π³ℏ / (G_F² m_μ⁵);   τ_τ = τ_μ · BR · e^(−10 S),  S = ½ ln(m_τ/m_μ)"],
      f"G_F = {N(GF)} MeV⁻², m_μ = {N(mmu,7)} MeV, m_τ = {N(mtau,6)} MeV, BR(τ→eνν) = {N(V('BR_tau_e'))} → S = {N(S_tm,5)}",
      [f"Muon lifetime from Fermi theory = {N(tau_mu_c*1,5)} s (measured {N(V('tau_mu_s'),7)} s: {N((tau_mu_c/V('tau_mu_s')-1)*100,2)}% off)",
       f"Tau lifetime = {N(tau_tau_c,4)} fs", rec_line("A6 tau lifetime", tau_tau_c, CL("A6"))],
      f"{N(V('tau_tau_fs'))} fs — {SRC('tau_tau_fs')}", f"within {N(V('tol_percent_1')*100,2)}%", j_rel(tau_tau_c, V("tau_tau_fs"), V("tol_percent_1")),
      "e^(−10 S) is the same as (m_μ/m_τ)⁵, the Fermi-theory mass scaling, so the exponent is not a fitted number.")
# anomaly sums (used by A7, B5, M16)
Yv = {k: fr_ for k, fr_ in V("anomaly_Y").items()}
from fractions import Fraction as Fr
Y = {k: Fr(a, b) for k, (a, b) in Yv.items()}
cond = {"SU(3)²U(1)": 2 * Y["YQ"] - Y["Yu"] - Y["Yd"], "SU(2)²U(1)": 3 * Y["YQ"] + Y["YL"],
        "gravity²U(1)": 6 * Y["YQ"] - 3 * Y["Yu"] - 3 * Y["Yd"] + 2 * Y["YL"] - Y["Ye"],
        "U(1)³": 6 * Y["YQ"] ** 3 - 3 * Y["Yu"] ** 3 - 3 * Y["Yd"] ** 3 + 2 * Y["YL"] ** 3 - Y["Ye"] ** 3}
anom_ok = all(v == 0 for v in cond.values())
cols = {n: n * (Fr(2, 3) - Fr(1, 3)) + (0 - 1) for n in V("colours_tested")}
only3 = [n for n, v in cols.items() if v == 0]
anom_lines = [f"{k} = {N(float(v),3)}" for k, v in cond.items()]
col_line = "Sum of electric charges in one generation for N colours: " + "; ".join(f"N = {I(n)} → {N(float(v),3)}" for n, v in cols.items())
check("A7", "A", "Quark charges and three colours (anomaly cancellation)", (R + ANCH["M16"], "record M16"),
      ["2Y_Q − Y_u − Y_d = 0;  3Y_Q + Y_L = 0;  6Y_Q − 3Y_u − 3Y_d + 2Y_L − Y_e = 0;  6Y_Q³ − 3Y_u³ − 3Y_d³ + 2Y_L³ − Y_e³ = 0",
       "charge sum per generation = N (Q_u + Q_d) + Q_e + Q_ν = 0"],
      "hypercharges Y_Q = " + fr(*V("anomaly_Y")["YQ"]) + ", Y_u = " + fr(*V("anomaly_Y")["Yu"]) + ", Y_d = " + fr(*V("anomaly_Y")["Yd"]) + ", Y_L = " + fr(*V("anomaly_Y")["YL"]) + ", Y_e = " + fr(*V("anomaly_Y")["Ye"]),
      anom_lines + [col_line, f"Only N = {I(only3[0]) if len(only3)==1 else '?'} gives zero."],
      "all four anomaly sums must be exactly zero (Geng and Marshak 1989; Minahan, Ramond and Warner 1990)", "exact rational arithmetic: zero means zero",
      "PASS" if anom_ok and len(only3) == 1 else "FAIL")
# ringdown
Msrc, zz, aa = V("ringdown_M_source"), V("ringdown_z"), V("ringdown_spin"); k = V("kerr_fit"); Msun_s = V("Msun_s")
Mrem = Msrc * (1 + zz)
Mw = k[0] - k[1] * (1 - aa) ** k[2]; Qf = k[3] + k[4] * (1 - aa) ** (-k[5])
w_ = Mw / (Mrem * Msun_s); f_ = w_ / (2 * math.pi); tau_ = 2 * Qf / w_ * 1e3
dt_echo = 4 * Mrem * Msun_s * (1 + 1 / math.sqrt(1 - aa * aa)) * math.log(Mrem * V("Msun_m") / V("planck_length_m"))
sf, st = V("ringdown_f_sigma"), V("ringdown_tau_sigma")
check("A8a", "A", "Black-hole ringdown tone and fade (GW150914)", (R + ANCH["K9"], "record K9"),
      ["Mω = a₀ − a₁ (1−χ)^a₂;  Q = b₀ + b₁ (1−χ)^(−b₂);  f = ω/2π;  τ = 2Q/ω"],
      f"remnant mass = {N(Msrc)} × (1 + z) = {N(Mrem,4)} solar masses (z = {N(zz)}), spin χ = {N(aa)}; fit constants {', '.join(N(v,5) for v in k)} (Berti, Cardoso and Will 2006)",
      [f"f = {N(f_,3)} Hz, τ = {N(tau_,3)} ms", rec_line("A8 frequency", f_, CL("A8f"), 3)],
      f"f = {N(V('ringdown_f_measured'))} ± {N(sf)} Hz, τ = {N(V('ringdown_tau_measured'))} ± {N(st)} ms — {SRC('ringdown_f_measured')}",
      "each within 1 standard deviation",
      "PASS" if j_sigma(f_, V("ringdown_f_measured"), sf) == "PASS" and j_sigma(tau_, V("ringdown_tau_measured"), st) == "PASS" else "FAIL")
check("A8b", "A", "No echoes after the ringdown", (R + ANCH["K9"], "record K9"),
      ["Δt = 4M(1 + 1/√(1−χ²)) ln(M/ℓ_P)   (delay if a reflecting surface sat one Planck length outside the horizon)"],
      f"M = {N(Mrem,4)} solar masses, χ = {N(aa)}, ℓ_P = {N(V('planck_length_m'))} m",
      [f"Echo delay that would appear = {N(dt_echo,3)} s. TRPG predicts there is no such surface, so no echoes."],
      "searches by LIGO/Virgo and by Westerweck and others (2018) found no significant echoes (this run does not repeat that search; the result is as the record cites it)",
      "a negative prediction is consistent with a null result", "PASS",
      "Consistent with the null searches. This is a prediction of absence, so it can still fail if echoes are found later.")
# neutrinos
dm21, dm31, dm32 = V("dm21"), V("dm31"), V("dm32_inverted")
def Qs(m, s): rt = np.array(s) * np.sqrt(m); return m.sum() / rt.sum() ** 2
def roots_of(f):
    grid = np.logspace(-6, 0, 4000); vals = [f(x) for x in grid]
    return [brentq(f, grid[i], grid[i + 1]) for i in range(len(grid) - 1) if np.sign(vals[i]) != np.sign(vals[i + 1])]
nrm = []
for sg in [(1, 1, 1), (-1, 1, 1)]:
    f = lambda m1, sg=sg: Qs(np.array([m1, math.sqrt(m1 ** 2 + dm21), math.sqrt(m1 ** 2 + dm31)]), sg) - Q23
    for m1 in roots_of(f): nrm.append((sg, m1, math.sqrt(m1 ** 2 + dm21), math.sqrt(m1 ** 2 + dm31)))
inv = []
for sg in [(1, 1, 1), (1, 1, -1), (1, -1, 1), (-1, 1, 1)]:
    def f(m3, sg=sg):
        m2 = math.sqrt(m3 ** 2 - dm32); m1 = math.sqrt(m2 ** 2 - dm21); return Qs(np.array([m1, m2, m3]), sg) - Q23
    for m3 in roots_of(f):
        m2 = math.sqrt(m3 ** 2 - dm32); inv.append((sg, math.sqrt(m2 ** 2 - dm21), m2, m3))
pos = [r for r in nrm if r[0] == (-1, 1, 1)]
pred = pos[0] if pos else None
sum_pred = sum(pred[1:]) if pred else float("nan")
lim = V("nu_sum_limit")
nrm_lines = [f"Normal ordering, square-root signs ({', '.join(('−' if s<0 else '+') for s in sg)}): m₁, m₂, m₃ = {N(a*1e3,3)}, {N(b*1e3,3)}, {N(c*1e3,3)} meV, sum = {N((a+b+c),4)} eV" for sg, a, b, c in nrm]
inv_lines = [f"Inverted ordering, signs ({', '.join(('−' if s<0 else '+') for s in sg)}): sum = {N(a+b+c,3)} eV" for sg, a, b, c in inv]
check("A9a", "A", "Neutrino masses from the lepton rule (normal ordering)", (R + ANCH["M2"], "record M2"),
      ["Q = (m₁ + m₂ + m₃)/(±√m₁ ± √m₂ ± √m₃)² = 2/3,  m₂² = m₁² + Δm²₂₁,  m₃² = m₁² + Δm²₃₁"],
      f"Δm²₂₁ = {N(dm21)} eV², Δm²₃₁ = {N(dm31)} eV² ({SRC('dm21')})",
      (nrm_lines + [rec_line("A9 neutrino sum", sum_pred, CL("A9")), rec_line("A9 lightest mass (meV)", pred[1] * 1e3, V("nu_light_stated_meV"), 3)]) if pred else nrm_lines,
      f"not yet measured; the record's stated prediction is {', '.join(N(v) for v in V('nu_masses_stated'))} meV", "the solution with one negative square root must reproduce the stated masses within " + N(V('tol_percent_5')*100,2) + "%",
      ("PASS" if pred and all(abs(p*1e3/s-1) <= V("tol_percent_5") for p, s in zip(pred[1:], V("nu_masses_stated"))) else "FAIL"),
      "This stage checks the arithmetic of the prediction, not nature. Nature's answer is in the next two stages.")
check("A9b", "A", "Does the predicted neutrino mass sum respect the cosmology limit?", (R + ANCH["M2"], "record M2"),
      ["Σm_ν (normal-ordering solution) < Σm_ν limit;   Σm_ν (inverted-ordering solutions) < limit?"],
      f"limit = {N(lim)} eV ({SRC('nu_sum_limit')}); predicted sum = {N(sum_pred,4)} eV",
      inv_lines if inv else ["no inverted-ordering solution found"],
      f"{N(lim)} eV upper limit — {SRC('nu_sum_limit')}", "predicted sum below the limit, and every inverted solution above it (ruling inverted ordering out)",
      "PASS" if sum_pred < lim and inv and all(sum(r[1:]) > lim for r in inv) else "FAIL")
num_min = math.sqrt(dm21) + math.sqrt(dm31)
check("A9d", "A", "How close the predicted sum is to the smallest sum normal ordering allows", (R + ANCH["M2"], "record M2"),
      ["Σm_min = √Δm²₂₁ + √Δm²₃₁   (lightest neutrino massless)"], f"Δm²₂₁ = {N(dm21)} eV², Δm²₃₁ = {N(dm31)} eV²",
      [f"Smallest allowed sum = {N(num_min,4)} eV; predicted sum = {N(sum_pred,4)} eV; limit = {N(lim)} eV", rec_line("A9 minimum sum", num_min, CL("numin"), 4)],
      "the record: the sum is only just above the smallest allowed and just under the limit", "predicted sum between the smallest allowed sum and the limit", j_range(sum_pred, num_min, lim),
      "The distinctive TRPG number is the lightest mass, not the sum.")
check("A9c", "A", "Neutrino ordering: normal, as predicted", (R + ANCH["M2"], "record M2"),
      ["prediction: normal ordering (m₁ < m₂ < m₃); inverted ordering would falsify it"],
      "no data yet", ["Awaiting the JUNO reactor experiment (running since 2025) and the sky surveys."],
      "not yet measured", "the JUNO measurement of the ordering", "NOT TESTABLE YET",
      "This is the prediction that could prove TRPG wrong. If JUNO finds inverted ordering, or a sum clearly away from the predicted one, this piece of TRPG fails.")

# ================================================================================================ SECTION B: master equation terms
SECTIONS["B"] = "The master equation, term by term"
# area law on the plain lattice
def lattice(Nn, m, crossing=False, amp=0.5):
    idx = lambda i, j: i * Nn + j; K = np.zeros((Nn * Nn, Nn * Nn))
    for i in range(Nn):
        for j in range(Nn):
            a = idx(i, j)
            for di, dj in ((1, 0), (0, 1)):
                ii, jj = i + di, j + dj
                if ii < Nn and jj < Nn:
                    xm, ym = (i + ii) / 2, (j + jj) / 2
                    kk = 1 + amp * math.cos(2 * math.pi * xm / 6) * math.cos(2 * math.pi * ym / 6) if crossing else 1.0
                    b = idx(ii, jj); K[a, b] = K[b, a] = -kk; K[a, a] += kk; K[b, b] += kk
            K[a, a] += m * m
    w, U = np.linalg.eigh(K); X = 0.5 * (U / np.sqrt(w)) @ U.T; P = 0.5 * (U * np.sqrt(w)) @ U.T
    return idx, X, P
def lattice_area(Nn, m):
    idx = lambda i, j: i * Nn + j; K = np.zeros((Nn * Nn, Nn * Nn))
    for i in range(Nn):
        for j in range(Nn):
            a = idx(i, j); K[a, a] = 4 + m * m
            for di, dj in ((1, 0), (0, 1)):
                ii, jj = i + di, j + dj
                if ii < Nn and jj < Nn: b = idx(ii, jj); K[a, b] = K[b, a] = -1
    w, U = np.linalg.eigh(K); X = 0.5 * (U / np.sqrt(w)) @ U.T; P = 0.5 * (U * np.sqrt(w)) @ U.T
    return idx, X, P
def ent(Nn, m, crossing, L, amp=0.5, area_script=False):
    idx, X, P = lattice_area(Nn, m) if area_script else lattice(Nn, m, crossing, amp); out = []
    for l in L:
        o = (Nn - l) // 2; A = [idx(i, j) for i in range(o, o + l) for j in range(o, o + l)]
        nu = np.sqrt(np.clip(np.linalg.eigvals(X[np.ix_(A, A)] @ P[np.ix_(A, A)]).real, 0.25 + 1e-12, None))
        out.append(sum((n + .5) * math.log(n + .5) - (n - .5) * math.log(n - .5) for n in nu if n > 0.5 + 1e-9))
    return np.array(out)
Ng, mg = int(V("grid_N")), V("grid_mass"); Ls = np.array(V("area_blocks"))
S_area = ent(Ng, mg, False, Ls, area_script=True)
cb = np.polyfit(4 * Ls, S_area, 1); cv = np.polyfit(Ls * Ls, S_area, 1)
rb = float(np.sqrt(np.mean((np.polyval(cb, 4 * Ls) - S_area) ** 2))); rv = float(np.sqrt(np.mean((np.polyval(cv, Ls * Ls) - S_area) ** 2)))
Gcoef = 1 / (4 * cb[0])
check("B1a", "B", "Gravity term: the area law holds in the wave medium", (R + ANCH["M16"], "record M16"),
      ["S(region) = c × (boundary length)   rather than   c × (area inside)",
       "G = a² / (4c)   (Bekenstein–Jacobson relation)"],
      f"{I(Ng)}×{I(Ng)} lattice of coupled oscillators, mass {N(mg)}, square blocks of side {', '.join(I(l) for l in Ls)}",
      [f"Entanglement per unit boundary: {', '.join(N(s/(4*l),3) for s, l in zip(S_area, Ls))}",
       f"Boundary-fit error = {N(rb,3)}; volume-fit error = {N(rv,3)}; ratio = {N(rv/rb,3)}", f"c = {N(cb[0],3)}, so G = {N(Gcoef,3)} a²",
       rec_line("B1 G coefficient", Gcoef, CL("B1"), 3)],
      "Bombelli, Koul, Lee and Sorkin (1986); Srednicki (1993): the area law is the known result for this kind of system",
      f"the boundary fit must beat the volume fit by at least a factor {N(V('area_ratio_min'))}", "PASS" if rv / rb >= V("area_ratio_min") else "FAIL")
check("B1b", "B", "Gravity term: Newton's G in real units", (R + ANCH["P1"], "record P1"),
      ["G = a² × (the coefficient above), with a = the medium's smallest spacing"],
      "a is not fixed by any measurement; the record takes it to be the Planck length",
      ["Not computed: the coefficient is in units of a², and a is not derived."], "G is measured, but TRPG does not yet predict a", "a derived value of a", "NOT TESTABLE YET")
Hs = V("H0") * 1e3 / V("Mpc_m"); rho_c = 3 * Hs ** 2 / (8 * math.pi * G)
uL = V("Omega_L_stated") * rho_c * V("c_exact") ** 2; uL_MeVfm3 = uL / V("MeV_J") * 1e-45
EL = (uL_MeVfm3 * V("hbarc_MeVfm") ** 3) ** 0.25 * 1e9
check("B2a", "B", "Cosmic baseline resistance −ρ_Λ", (R + ANCH["M8"], "record M8"),
      ["ρ_c = 3H₀²/(8πG);   u_Λ = Ω_Λ ρ_c c²;   E = (u_Λ (ℏc)³)^(1/4)"],
      f"H₀ = {N(V('H0'))} km/s/Mpc, Ω_Λ = {N(V('Omega_L_stated'))} ({SRC('H0')})",
      [f"u_Λ = {N(uL_MeVfm3,3)} MeV/fm³; energy scale = {N(EL,3)} meV", rec_line("B2 density", uL_MeVfm3, V("de_density_stated"), 3), rec_line("B2 scale", EL, V("de_scale_meV_stated"), 3)],
      "taken from the measured dark-energy density: an input, not a TRPG prediction", "not applicable", "INPUT",
      "TRPG does not yet derive this value.")
nu_lo, nu_hi = min(V("nu_masses_stated")), max(V("nu_masses_stated"))
check("B2b", "B", "Dark-energy scale versus the neutrino masses", (R + ANCH["M8"], "record M8"),
      ["m_ν,min ≤ E_Λ ≤ m_ν,max"],
      f"E_Λ = {N(EL,3)} meV; predicted neutrino masses {N(nu_lo)} to {N(nu_hi)} meV",
      [f"E_Λ lies {'inside' if nu_lo <= EL <= nu_hi else 'outside'} the predicted range."], "known puzzle: Fardon, Nelson and Weiner (2004)", "inside the range",
      j_range(EL, nu_lo, nu_hi), "A coincidence of scale, not a derivation.")
hbc = V("hbarc_MeVfm"); B_list = V("bag_B_proton")
gl_list = [V("glueball_lo"), V("glueball_hi")]
rows_lv = []
for Bv in B_list:
    for mgl in gl_list:
        Bn = Bv * hbc ** 3; lam = mgl ** 4 / (16 * Bn); v_ = math.sqrt(mgl ** 2 / (2 * lam)); rows_lv.append((Bv, mgl, lam, v_, hbc / mgl))
lam_lo = min(r[2] for r in rows_lv); lam_hi = max(r[2] for r in rows_lv); v_lo = min(r[3] for r in rows_lv); v_hi = max(r[3] for r in rows_lv)
wall = rows_lv[0][4]; core_hi = V("core_radius_hi")
cen_rows = [r for r in rows_lv if r[0] == B_list[1]]; lam_c_lo = min(r[2] for r in cen_rows); lam_c_hi = max(r[2] for r in cen_rows); v_c = [r[3] for r in cen_rows]
check("B3a", "B", "Containment field P: stiffness and depth", (R + ANCH["M1"], "record M1"),
      ["V(P) = (λ/4) v⁴ (P² − 1)²;   depth V(0) − V(1) = B;   glueball mass m² = 2λv²;   so λ = m⁴/(16B),  v = m/√(2λ)"],
      "B = " + ", ".join(N(b) for b in B_list) + " MeV/fm³ (set by the proton, stage B9); glueball = " + ", ".join(N(m) for m in gl_list) + " MeV (Morningstar and Peardon 1999; Chen and others 2006)",
      [f"λ from {N(lam_lo,3)} to {N(lam_hi,3)}; v from {N(v_lo,3)} to {N(v_hi,3)} MeV", f"Wall thickness ≈ ℏc/m = {N(wall,3)} fm",
       f"At the central B = {N(B_list[1])}: λ from {N(lam_c_lo,3)} to {N(lam_c_hi,3)}; v from {N(min(v_c),3)} to {N(max(v_c),3)} MeV",
       rec_line("B3 lambda_low", lam_c_lo, CL("B3lo"), 3), rec_line("B3 lambda_high", lam_c_hi, CL("B3hi"), 3), rec_line("B3 v", sum(v_c) / len(v_c), CL("B3v"), 3)],
      "no measurement of λ or v exists; this stage checks that the wall is thin enough for the bag-model results to carry over",
      "wall thickness must be well below the proton core radius " + N(core_hi) + " fm", "PASS" if wall < 0.5 * core_hi else "FAIL",
      "The quartic shape (P² − 1)² is the simplest choice, not a derived one (record M1).")
eps = (V("gluon_vev_eps_prefactor_num") / V("gluon_vev_eps_prefactor_den")) * V("svz_G2") * 1e12 / hbc ** 3
ratio_svz = eps / B_list[1]
check("B3b", "B", "Containment depth versus the standard QCD vacuum energy", (R + ANCH["M1"], "record M1"),
      ["ε_SVZ = (9/32) ⟨α_s G²/π⟩;   compare B"],
      f"⟨α_s G²/π⟩ = {N(V('svz_G2'))} GeV⁴ ({SRC('svz_G2')}); B = {N(B_list[1])} MeV/fm³",
      [f"ε_SVZ = {N(eps,3)} MeV/fm³, which is {N(ratio_svz,3)} times B"],
      f"the standard QCD estimate is {N(eps,3)} MeV/fm³", f"B within a factor {N(V('svz_factor'))} of the QCD estimate",
      "PASS" if ratio_svz <= V("svz_factor") else "FAIL", "The record calls this a known, unresolved gap in mainstream physics. It is shown here as a FAIL because it is one.")
rng_ = np.random.default_rng(int(V("random_seed")))
mvals = rng_.uniform(0.1, 100, 6)
scale_f = V("scale_factors"); S_ok = True; cosh_ok = True
for lam_ in scale_f:
    for i_ in range(5):
        a1, a2 = mvals[i_], mvals[i_ + 1]
        S1 = 0.5 * math.log(a1 / a2); S2 = 0.5 * math.log(lam_ * a1 / (lam_ * a2))
        if abs(S1 - S2) > 1e-12: S_ok = False
        if abs((a1 + a2) / math.sqrt(a1 * a2) - 2 * math.cosh(S1)) > 1e-9: cosh_ok = False
check("B4", "B", "Scale layer σ: the distance S is the only one that gives no particle a privileged unit", (R + ANCH["M3"], "record M3"),
      ["S_ij = ½ ln(m_i/m_j);  S(λm_i, λm_j) = S(m_i, m_j);  M_ij = (m_i + m_j)/√(m_i m_j) = 2 cosh S_ij"],
      f"random positive masses; rescaling factors λ = {', '.join(N(v) for v in scale_f)}",
      [f"S unchanged by rescaling in every trial: {'yes' if S_ok else 'no'}", f"M_ij = 2 cosh S_ij in every trial: {'yes' if cosh_ok else 'no'}"],
      "exact identities", f"agreement to {N(V('tol_identity'))} relative", "PASS" if S_ok and cosh_ok else "FAIL",
      "This checks the algebra. That S is forced by the symmetry is the record's argument (M3).")
check("B5", "B", "Forces: charges and colours", (R + ANCH["M16"], "record M16"),
      ["−¼ F² for U(1) × SU(2) × SU(3);  charges forced by anomaly cancellation"],
      "see stage A7 for the sums", [f"Anomaly sums all zero: {'yes' if anom_ok else 'no'}.", "The force groups themselves are chosen, not derived: the record says so (T2 is PARTIAL)."],
      "anomaly cancellation (Geng and Marshak 1989; Minahan, Ramond and Warner 1990)", "exact zero", "PARTIAL",
      "Charges and colour count are forced; the groups U(1), SU(2), SU(3) are not.")
check("B6", "B", "Matter term ψ̄ iγDψ", (R + ANCH["P1"], "record P1"),
      ["ψ̄ iγ^μ D_μ ψ"], "Dirac (1928)", ["Taken unchanged from established physics; TRPG adds no number here."], "Dirac (1928)", "not applicable", "INHERITED")
# lepton matrix
r_ = 1 / math.sqrt(2); d_ = D29; M0 = (math.sqrt(me) + math.sqrt(mmu) + math.sqrt(mtau)) / 3
w_c = np.exp(1j * d_); Sm = M0 * np.array([[1, r_ * w_c, r_ * np.conj(w_c)], [r_ * np.conj(w_c), 1, r_ * w_c], [r_ * w_c, r_ * np.conj(w_c), 1]])
ev = np.sort(np.linalg.eigvalsh(Sm)) ** 2
vv = np.sqrt([me, mmu, mtau]); cc = vv.mean() * np.ones(3); pp = vv - cc
errs = [abs(ev[0] / me - 1), abs(ev[1] / mmu - 1), abs(ev[2] / mtau - 1)]
check("B7a", "B", "Mass couplings: the lepton matrix gives the three masses", (R + ANCH["M2"], "record M2"),
      ["S_ℓ = M₀ × circulant(1, e^{iδ}/√2, e^{−iδ}/√2);  masses = (eigenvalues of S_ℓ)²"],
      f"δ = {fr(V('brannen_delta_num'), V('brannen_delta_den'))} = {N(D29,6)}; M₀ = (√m_e + √m_μ + √m_τ)/3 = {N(M0,6)} MeV½",
      [f"Eigenvalues squared: {N(ev[0],6)}, {N(ev[1],7)}, {N(ev[2],6)} MeV"],
      f"{N(me,6)}, {N(mmu,7)}, {N(mtau,6)} MeV — {SRC('m_e')}", f"each within {N(V('tol_percent_1')*100,2)}%",
      "PASS" if all(e <= V("tol_percent_1") for e in errs) else "FAIL",
      "M₀ is computed from the measured masses here, so this stage shows the matrix form is consistent with the three masses, not that it predicts them from nothing.")
check("B7b", "B", "Self part equals relational part (Koide's balance)", (R + ANCH["M2"], "record M2"),
      ["|common part|² = |relational part|² of the square-root-mass vector"],
      f"common part² = {N(cc@cc,6)}, relational part² = {N(pp@pp,6)} MeV",
      [f"Ratio = {N((pp@pp)/(cc@cc),6)}; common part² versus the proton mass: {N((cc@cc/V('m_proton')-1)*100,3)}%"],
      "Koide (1981); Foot (1994)", f"ratio equal to 1 within {N(V('tol_lepton_mass')*100,2)}%", j_rel((pp@pp) / (cc@cc), 1.0, V("tol_lepton_mass")))
# ================================================================================================ SECTION C: chapter K
SECTIONS["C"] = "Chapter K: the fixes, one by one"
gm_mmax, gm_r14, gm_p150 = gemini_sequence()
check("K1a", "C", "Gemini's own neutron-star equation solved properly", (R + ANCH["K1"], "record K1"),
      ["dp/dr = −(ε + p)(m + 4πr³p)/(r(r − 2m)),  dm/dr = 4πr²ε  with Gemini's low-density formula"],
      f"B = {N(V('bag_B_gemini'))} MeV/fm³, low-density term p = {N(V('gemini_p_scale'))} ε (ε/4B)^{N(V('gemini_p_power'))}",
      [f"R({N(STDM)}) = {N(gm_r14,3)} km; maximum mass {N(gm_mmax,3)} solar masses", f"Pressure at ε = {N(V("nuclear_density_probe"))} MeV/fm³ is {N(gm_p150,3)} MeV/fm³ (real nuclear matter is about {N(V('nuclear_pressure_range')[0])} to {N(V('nuclear_pressure_range')[1])})", rec_line("K1a R1.4", gm_r14, CL("K1aR"), 3)],
      f"{N(V('ns_radius_lo'))} to {N(V('ns_radius_hi'))} km — {SRC('ns_radius_lo')}", "inside the measured range", j_range(gm_r14, V("ns_radius_lo"), V("ns_radius_hi")),
      "This is the version that was wrong. It is shown so the failure is on the page.")
Es, Ps = sly_tables()
sly_pf = lambda e: float(np.interp(e, Es, Ps)); sly_inv = lambda p: float(np.interp(p, Ps, Es))
sl = np.array([(ec,) + tov(sly_pf, sly_inv, ec, psurf=1e-12) for ec in np.linspace(V("sly_scan")[0], V("sly_scan")[1], 57)])
i_ = sl[:, 2].argmax(); jj = int(np.argmin(abs(sl[:i_, 2] - STDM)))
sly_m, sly_r = float(sl[i_, 2]), float(sl[jj, 1])
check("K1b", "C", "Solver check on a standard real equation of state (SLy)", (R + ANCH["K1"], "record K1"),
      ["same TOV equations, SLy piecewise-polytrope equation of state"],
      "SLy parameters (Read, Lackey, Owen and Friedman 2009), central density scanned over " + N(V("sly_scan")[0]) + " to " + N(V("sly_scan")[1]) + " MeV/fm³",
      [f"Maximum mass = {N(sly_m,3)} solar masses; R({N(STDM)}) = {N(sly_r,3)} km", rec_line("K1b SLy mass", sly_m, V("sly_published_mmax"), 3)],
      f"{N(V('sly_published_mmax'))} solar masses and {N(V('sly_published_r14'))} km — {SRC('sly_published_mmax')}", f"each within {N(V('tol_percent_5')*100,2)}%",
      "PASS" if j_rel(sly_m, V("sly_published_mmax"), V("tol_percent_5")) == "PASS" and j_rel(sly_r, V("sly_published_r14"), V("tol_percent_5")) == "PASS" else "FAIL",
      "If this solver were wrong, none of the other star results could be trusted.")
mmw, r14w, _ = bag_sequence(bB, V("bag_wall_gemini"))
check("K1c", "C", "Containment-only star, and with Gemini's wall term", (R + ANCH["K1"], "record K1"),
      ["p = (ε − 4B)/3 + wall term; same TOV equations"],
      f"B = {N(bB)} MeV/fm³; wall strength = {N(V('bag_wall_gemini'))}",
      [f"Containment only: M_max = {N(mm,3)}, R({N(STDM)}) = {N(r14,3)} km", f"With wall term: M_max = {N(mmw,3)}, R({N(STDM)}) = {N(r14w,3)} km", rec_line("K1c wall mass", mmw, CL("K1cMw"), 3)],
      f"{N(V('ns_mass_measured'))} ± {N(V('ns_mass_sigma'))} solar masses; {N(V('ns_radius_lo'))} to {N(V('ns_radius_hi'))} km", "heaviest star within 1 standard deviation or above it, and radius in range",
      "PASS" if mm >= V("ns_mass_measured") - V("ns_mass_sigma") and mmw >= V("ns_mass_measured") - V("ns_mass_sigma") and V("ns_radius_lo") <= r14 <= V("ns_radius_hi") and V("ns_radius_lo") <= r14w <= V("ns_radius_hi") else "FAIL")
fj4 = (bB * hbc ** 3) ** 0.25; fj_min4 = (V("bag_B_fj_min") * hbc ** 3) ** 0.25
check("K1d", "C", "The stability rule: is Gemini's B allowed?", (R + ANCH["K1"], "record K1"),
      ["B ≥ B_min  (Farhi and Jaffe 1984), i.e. B¹ᐟ⁴ ≥ about 145 MeV"],
      f"B = {N(bB)} MeV/fm³ → B^(1/4) = {N(fj4,4)} MeV; the bound B_min = {N(V('bag_B_fj_min'))} MeV/fm³ → {N(fj_min4,4)} MeV",
      [f"B = {N(bB)} is {'above' if bB >= V('bag_B_fj_min') else 'below'} the bound."],
      f"B_min = {N(V('bag_B_fj_min'))} MeV/fm³ — {SRC('bag_B_fj_min')}", "B at or above B_min", "PASS" if bB >= V("bag_B_fj_min") else "FAIL",
      "The record says this too: this B is slightly too small, and the fix is a larger B.")
# K2 tidal
cs2_arr = np.gradient(Ps, Es)
dpde_s = lambda e: float(np.interp(e, Es, cs2_arr))
sl2 = np.array([tidal_star(sly_pf, sly_inv, dpde_s, ec, 0.0, psurf=1e-9) for ec in np.linspace(V("sly_scan")[0], V("sly_scan")[1], 120)])
ii = sl2[:, 0].argmax(); rr_ = sl2[:ii + 1]
sly_lam = float(np.exp(np.interp(STDM, rr_[:, 0], np.log(rr_[:, 3]))))
check("K2a", "C", "Solver check for the tidal number Λ (SLy)", (R + ANCH["K2"], "record K2"),
      ["k₂ from y(R) and C = M/R (Hinderer 2008);  Λ = (2/3) k₂ / C⁵"],
      f"SLy equation of state; star of {N(STDM)} solar masses",
      [f"Λ({N(STDM)}) = {N(sly_lam,3)}", rec_line("K2a SLy Lambda", sly_lam, CL("K2a"), 3)],
      f"about {N(V('sly_published_lambda'))} — {SRC('sly_published_lambda')}", f"within {N(V('tol_percent_5')*100,2)}%", j_rel(sly_lam, V("sly_published_lambda"), V("tol_percent_5")))
tid_lines = [ln for ln in PUB["neutron_star_tidal_fix"][2].splitlines() if "Lambda(" + str(STDM) + ")" in ln]
for n_, ln in enumerate(tid_lines):
    Bv, wl = re.findall(r"B=([\d.]+) wall=([\d.]+)", ln)[0]; Mx = float(re.search(r"Mmax=([\d.]+)", ln).group(1)); Lm = float(re.search(r"Lambda\(" + re.escape(str(STDM)) + r"\)=(\d+)", ln).group(1))
    ok2 = Mx >= mt and Lm <= V("lambda_limit")
    check(f"K2-{n_ + 1}", "C", f"Tidal number from the published script, run {I(n_ + 1)}", (R + ANCH["K2"], "record K2"),
          ["k\u2082 and \u039b from the TOV and Love-number equations (Hinderer 2008), self-bound surface correction"],
          f"B = {N(float(Bv))} MeV/fm\u00b3, wall strength {N(float(wl))} (printed by neutron_star_tidal_fix)",
          [f"M_max = {N(Mx,3)} solar masses; \u039b({N(STDM)}) = {N(Lm,3)}"], f"M_max at least {N(mt)}; \u039b at most {N(V('lambda_limit'))} \u2014 {SRC('lambda_limit')}",
          "both tests", "PASS" if ok2 else "FAIL",
          "These runs fit with two adjustable numbers (B and the wall strength): superseded in the record by the no-knobs result (stage A3)." if float(wl) > 0 else "")
# K3
mp_ = V("m_proton")
def B_of_R(Rf): return 3 * mp_ / (16 * math.pi * Rf ** 3)
B_full = B_of_R(V("proton_radius_fm")); Bc_lo = B_of_R(V("core_radius_hi")); Bc_hi = B_of_R(V("core_radius_lo"))
def R_of_B(Bv): return (3 * mp_ / (16 * math.pi * Bv)) ** (1 / 3)
R_in_lo, R_in_hi = R_of_B(V("bag_B_window_hi")), R_of_B(V("bag_B_window_lo"))
rms_f = math.sqrt(3 / 5)
check("K3a", "C", "Proton containment balance: B from the proton's whole radius", (R + ANCH["K3"], "record K3"),
      ["M = 4 B V,  V = (4/3)πR³   →   B = 3M/(16πR³)"],
      f"M = {N(mp_,6)} MeV; R = {N(V('proton_radius_fm'))} fm (the whole measured radius)",
      [f"B = {N(B_full,3)} MeV/fm³", rec_line("K3 B whole radius", B_full, CL("K3B"), 3),
       f"For comparison the stability bound is {N(V('bag_B_fj_min'))} MeV/fm³, so this B is {'above' if B_full >= V('bag_B_fj_min') else 'below'} it."],
      f"the record states about {N(CL('K3B'))} MeV/fm³, 'too small' (below the stability bound {N(V('bag_B_fj_min'))})", f"within {N(V('tol_percent_5')*100,2)}% of the record's figure",
      j_rel(B_full, CL("K3B"), V("tol_percent_5")),
      "The record's own formula does not give its stated number, and the result is the opposite side of the stability bound. This is flagged for the authors to reconcile.")
core_rms_lo, core_rms_hi = rms_f * R_in_lo, rms_f * R_in_hi
dev = max(abs(core_rms_lo / V("core_radius_lo") - 1), abs(core_rms_hi / V("core_radius_hi") - 1))
check("K3c", "C", "The quark-core size that gives a B inside the neutron-star window", (R + ANCH["K3"], "record K3"),
      ["R = (3M / (16πB))^(1/3);   rms radius of a uniform sphere = √(3/5) R"],
      f"M = {N(mp_,6)} MeV; B from {N(V('bag_B_window_lo'))} to {N(V('bag_B_window_hi'))} MeV/fm³ (the window)",
      [f"R from {N(R_in_lo,3)} to {N(R_in_hi,3)} fm; as an rms radius {N(core_rms_lo,3)} to {N(core_rms_hi,3)} fm",
       f"The record states a core of {N(V('core_radius_lo'))} to {N(V('core_radius_hi'))} fm (rms)."],
      f"{N(V('core_radius_lo'))} to {N(V('core_radius_hi'))} fm — {SRC('core_radius_lo')}", f"both ends within {N(V('tol_percent_5')*100,2)}% of the record's range",
      "PASS" if dev <= V("tol_percent_5") else ("PARTIAL" if dev <= PF * V("tol_percent_5") else "FAIL"),
      "The rms conversion (uniform sphere) is this page's reading of the record's 'rms'. The record's M1 text uses the bag radius itself, which matches the R computed here.")
check("K3b", "C", "Gemini's centroid and its 'hyperfine' number", (R + ANCH["K3"], "record K3"),
      ["centroid = (M_N + M_Δ)/2;  C_ss = gap / 6;  predicted gap = 6 C_ss"],
      f"M_N = {N(V('nucleon_mass_avg'))} MeV, M_Δ = {N(V('delta_mass'))} MeV; measured gap = {N(V('np_delta_gap'))} MeV; C_ss = {N(V('gemini_C_ss'))} MeV",
      [f"Centroid = {N((V('nucleon_mass_avg')+V('delta_mass'))/2,5)} MeV; gap / 6 = {N(V('np_delta_gap')/6,4)} MeV; 6 × C_ss = {N(6*V('gemini_C_ss'),5)} MeV"],
      f"measured gap {N(V('np_delta_gap'))} MeV", "a prediction must not contain its own answer",
      "FAIL" if abs(6 * V("gemini_C_ss") - V("np_delta_gap")) / V("np_delta_gap") < V("tol_percent_1") else "PASS",
      "Circular: C_ss was typed in as the measured gap divided by 6, so the 'prediction' returns the input.")
# K4 other mixings
pairs = [("V_ub", "q_u", "q_t", "ckm_Vub"), ("V_cb", "q_c", "q_t", "ckm_Vcb"), ("V_td", "q_d", "q_b", "ckm_Vtd"), ("V_ts", "q_s", "q_b", "ckm_Vts")]
for nm, a, b, ck in pairs:
    val = math.sqrt(V(a) / V(b)); meas = V(ck)
    check("K4-" + nm[2:], "C", f"The same √(m_i/m_j) pattern for {nm}", (R + ANCH["K4"], "record K4"),
          [f"{nm} ≈ √(m_{a[2]} / m_{b[2]})"], f"m_{a[2]} = {N(V(a))} MeV, m_{b[2]} = {N(V(b))} MeV",
          [f"√(m_{a[2]}/m_{b[2]}) = {N(val,3)}; ratio to measured = {N(val/meas,3)}"],
          f"{nm} = {N(meas)} — {SRC(ck)}", f"within {N(V('tol_percent_5')*100,2)}%", j_rel(val, meas, V("tol_percent_5")),
          "The record names V_ub explicitly and says V_cb, V_td and V_ts fail; the mass pairs for those three are the ones the same pattern points to, chosen here." if nm != "V_ub" else "The record: roughly works for V_ub, a little low.")
# K5 / K10 counts
if XLSX.exists():
    parts = read_xlsx(XLSX, "Particles"); pairs_x = read_xlsx(XLSX, "All_Pairs")
    cnt = [len(parts), len(pairs_x), sum(1 for r in pairs_x if (r.get("charge_product") or 0) != 0), sum(1 for r in pairs_x if r.get("strong_link") == 1),
           sum(1 for r in pairs_x if (r.get("W_transition_amplitude") or 0) > 0)]
    st5 = V("xlsx_stated_counts")
    check("K5", "C", "Every-particle-to-every-particle table: the counts", (R + ANCH["K5"], "record K5"),
          ["count rows of the downloadable workbook"], "TRPG_Every_Particle_to_Every_Particle.xlsx, sheets Particles and All_Pairs (read by this script)",
          [f"Particles = {I(cnt[0])}; pairs = {I(cnt[1])}; linked by the electric force = {I(cnt[2])}; by the strong force = {I(cnt[3])}; by weak transitions = {I(cnt[4])}"],
          "the record states " + ", ".join(I(v) for v in st5) + " (in the same order)", "exact counts", "PASS" if cnt == st5 else "FAIL")
    pa = read_xlsx(XLSX, "Particle_vs_Antiparticle")
    fates = [r.get("fate") or "" for r in pa]
    fc = {"opposites": sum(int(r.get("direct_opposites") or 0) for r in pa), "annihilate": sum(1 for f_ in fates if f_.startswith("can annihilate -> neutral")),
          "pions": sum(1 for f_ in fates if "pions" in f_), "W": sum(1 for f_ in fates if "fuse" in f_),
          "flavour": sum(1 for f_ in fates if f_.startswith("flavour-changing")), "scatter": sum(1 for f_ in fates if f_.startswith("no annihilation"))}
    sf5 = V("fate_stated")
    check("K10a", "C", "Particle-versus-antiparticle table: the counts (corrected rules)", (R + ANCH["K10"], "record K10"),
          ["count rows of the downloadable workbook"], "sheet Particle_vs_Antiparticle of the same workbook (read by this script)",
          [f"Direct opposites = {I(fc['opposites'])}; annihilate into neutral carriers = {I(fc['annihilate'])}; into pions = {I(fc['pions'])}; fuse into a W = {I(fc['W'])}; flavour-changing = {I(fc['flavour'])}; scatter only = {I(fc['scatter'])}"],
          "the record states " + ", ".join(f"{k} {I(v)}" for k, v in sf5.items()), "exact counts", "PASS" if fc == sf5 else "FAIL")
    # M15 network (force carriers included)
    P_ = {r["name"]: r for r in parts}
    def graph(use_higgs):
        E = {}
        for r in pairs_x:
            linked = (r.get("charge_product") or 0) != 0 or r.get("strong_link") == 1 or (r.get("W_transition_amplitude") or 0) > 0 or (use_higgs and (r.get("Higgs_coupling_product") or 0) > 0)
            if linked: E.setdefault(r["particle_i"], set()).add(r["particle_j"]); E.setdefault(r["particle_j"], set()).add(r["particle_i"])
        return E
    names = sorted(set(r["particle_i"] for r in pairs_x) | set(r["particle_j"] for r in pairs_x))
    E = graph(True)
    def link(a, b): E.setdefault(a, set()).add(b); E.setdefault(b, set()).add(a)
    for n_, r in P_.items():
        if n_ in ("photon", "gluon", "W", "Z", "Higgs"): continue
        if (r.get("charge") or 0) != 0: link("photon", n_)
        if r.get("colored"): link("gluon", n_)
        if (r.get("T3") or 0) != 0 or (r.get("charge") or 0) != 0: link("Z", n_); link("W", n_)
        if (r.get("mass_MeV") or 0) > 0: link("Higgs", n_)
    link("W", "photon"); link("W", "Z"); link("Higgs", "W"); link("Higgs", "Z")
    dist = []; unreached = 0
    for a in names:
        seen = {a: 0}; fr_ = [a]
        while fr_:
            nx = []
            for u in fr_:
                for v in E.get(u, ()):
                    if v not in seen: seen[v] = seen[u] + 1; nx.append(v)
            fr_ = nx
        for b in names:
            if b > a:
                if b in seen: dist.append(seen[b])
                else: unreached += 1
    nstat = V("network_stated")
    check("M15b", "D", "Degrees of separation in the particle table, with force carriers", (R + ANCH["M15"], "record M15"),
          ["shortest path over links (any direct interaction); average over all pairs"], "the particle workbook plus photon, gluon, W, Z and Higgs links",
          [f"Longest route = {I(max(dist))} steps; average = {N(sum(dist)/len(dist),3)}; pairs never connected = {I(unreached)}"],
          f"the record states longest {I(nstat[0])}, average {N(nstat[1])}", "same longest route and average within 1%",
          "PASS" if max(dist) == nstat[0] and abs(sum(dist) / len(dist) / nstat[1] - 1) <= V("tol_percent_1") and unreached == 0 else "FAIL")
# K6
Hs2 = V("H0") * 1e3 / V("Mpc_m"); rho_c2 = 3 * Hs2 ** 2 / (8 * math.pi * G); Rh = V("c_exact") / Hs2
Mu = (4 * math.pi / 3) * Rh ** 3 * rho_c2; mpk = V("m_proton_kg")
n_part = Mu / mpk
lam_c = V("hbarc_MeVfm") * 1e-15 / mp_
hbar_SI = V("hbarc_MeVfm") * 1e-15 * V("MeV_J") / V("c_exact")
alpha_g = G * mpk ** 2 / (hbar_SI * V("c_exact"))
wein = ((hbar_SI ** 2 * Hs2) / (G * V("c_exact"))) ** (1 / V("weinberg_relation_exponent_den"))
wein_MeV = wein * V("c_exact") ** 2 / V("MeV_J")
check("K6", "C", "Large-number coincidences in the proton's scale view", (R + ANCH["K6"], "record K6"),
      ["N = M_universe / m_p;  M_universe = (4π/3)(c/H₀)³ ρ_c;  (c/H₀)/ƒλ_p;  G m_p²/(ℏc);  m_W = (ℏ²H₀/(Gc))^(1/3)"],
      f"H₀ = {N(V('H0'))} km/s/Mpc, m_p = {N(mpk)} kg",
      [f"Particles in the universe's mass = {N(n_part,3)}", f"Universe across in proton Compton lengths = {N(Rh/(lam_c),3)}", f"Gravity weaker than the quantum scale by {N(1/alpha_g,3)}", f"Weinberg relation gives {N(wein_MeV,3)} MeV",
       rec_line("K6 N particles", n_part, CL("K6N"), 3)],
      f"Eddington's number is about {N(V('eddington_stated'))} (Eddington 1930s, Dirac 1937, Weinberg 1972)", "same order of magnitude as the stated number", "PARTIAL" if abs(math.log10(n_part / V("eddington_stated"))) < 1 else "FAIL",
      "PARTIAL because it is a coincidence, not a derivation, and the record notes every particle's view gives some huge number.")
# K7
rng2 = np.random.default_rng(int(V("random_seed"))); Nr = int(V("random_triples")); lo_, hi_ = V("random_range_log10")
mr = 10 ** rng2.uniform(lo_, hi_, (Nr, 3)); qr = mr.sum(1) / np.sqrt(mr).sum(1) ** 2
tolq = abs(Qlep - Q23); chance = float(np.mean(abs(qr - Q23) <= tolq))
check("K7a", "C", "Koide's relation on the charged leptons, and how special it is", (R + ANCH["K7"], "record K7"),
      ["Q = (m_e + m_μ + m_τ)/(√m_e + √m_μ + √m_τ)² = 2/3"],
      f"m_e, m_μ, m_τ as in stage A4; random test: {I(Nr)} log-uniform mass triples between {N(10**lo_)} and {N(10**hi_)} MeV",
      [f"Q = {N(Qlep,7)} (2/3 = {N(Q23,7)})", f"Chance of a random triple matching as closely = {N(chance,2)}", rec_line("K7 random chance", chance, V("random_chance_stated"), 2)],
      f"2/3 — Koide (1981)", f"Q within {N(V('tol_lepton_mass')*100,2)}% of 2/3", j_rel(Qlep, Q23, V("tol_lepton_mass")))
qk = {"u": V("q_u"), "d": V("q_d"), "s": V("q_s"), "c": V("q_c"), "b": V("q_b"), "t": V("q_t")}
for trip, nm in [("uct", "up-type (u, c, t)"), ("dsb", "down-type (d, s, b)"), ("cbt", "heavy trio (c, b, t)")]:
    qv = Qk(*[qk[x] for x in trip])
    check("K7-" + trip, "C", f"Koide's relation for quarks: {nm}", (R + ANCH["K7"], "record K7"),
          ["Q = Σm / (Σ√m)² = 2/3 ?"], "masses: " + ", ".join(f"{x} = {N(qk[x])} MeV" for x in trip),
          [f"Q = {N(qv,4)}"], "2/3", f"within {N(V('tol_percent_1')*100,2)}%", j_rel(qv, Q23, V("tol_percent_1")),
          "Quarks do not follow it as families; only the heavy trio comes near (record K7)." if trip != "cbt" else "")
check("K7-knobs", "C", "Gemini's mass formula: knobs versus results", (R + ANCH["K7"], "record K7"),
      ["knobs must be fewer than the masses they predict"], "the record's list: " + "; ".join(V("gemini_knob_count_items")),
      [f"Knobs listed = {I(len(V('gemini_knob_count_items')))}; masses predicted = {I(V('gemini_masses_count'))}"],
      "test T5 (fewer knobs than results)", "knobs fewer than results", "PASS" if len(V("gemini_knob_count_items")) < V("gemini_masses_count") else "FAIL")
check("K7-proton", "C", "Brannen scale versus the proton mass", (R + ANCH["K7"], "record K7"),
      ["3 M² ≈ m_p ?"], f"3M² = {N(3*Mt*Mt,5)} MeV; m_p = {N(mp_,6)} MeV",
      [f"Difference = {N((3*Mt*Mt/mp_-1)*100,2)}%"], f"{N(mp_,6)} MeV — {SRC('m_proton')}", f"within {N(V('proton_scale_hint')*100,2)}%",
      "PARTIAL" if abs(3 * Mt * Mt / mp_ - 1) <= V("proton_scale_hint") * 1.1 else "FAIL", "A hint only: unexplained, as the record says.")
# K8 Gemini claims
cent = V("gemini_centroid"); Css = V("gemini_C_ss"); MN, MD = cent - 3 * Css, cent + 3 * Css; mpi = V("pion_mass")
pp_ = math.sqrt((MD ** 2 - (MN + mpi) ** 2) * (MD ** 2 - (MN - mpi) ** 2)) / (2 * MD); EN = math.sqrt(MN ** 2 + pp_ ** 2)
Gam = lambda f: (f ** 2 / (12 * math.pi * mpi ** 2)) * ((EN + MN) / (2 * MD)) * pp_ ** 3
gw = Gam(V("gemini_f")); f_need = V("gemini_f") * math.sqrt(V("delta_width_measured") / gw)
check("K8a", "C", "Gemini's Δ decay width", (R + ANCH["K8"], "record K8"),
      ["Γ = f² p³ (E_N + M_N) / (12π m_π² · 2M_Δ)"],
      f"f = {N(V('gemini_f'))} (typed by Gemini), M_N = {N(MN,5)}, M_Δ = {N(MD,5)}, m_π = {N(mpi)} MeV",
      [f"Width = {N(gw,3)} MeV; Gemini claimed {N(V('gemini_width_claim'))} MeV; the f needed for the measured value is {N(f_need,3)}"],
      f"{N(V('delta_width_measured'))} MeV — {SRC('delta_width_measured')}", f"within {N(V('tol_percent_5')*100,2)}%", j_rel(gw, V("delta_width_measured"), V("tol_percent_5")),
      "The claimed number does not come out of Gemini's own formula.")
Z2 = 1 / V("gemini_Z2_denominator"); gA_f = Z2 * (V("gemini_gA_bare_sq") - (2 / 27) * V("gemini_gA_delta") * V("gemini_gA_bare_sq"))
check("K8b", "C", "Gemini's axial coupling g_A", (R + ANCH["K8"], "record K8"),
      ["g_A = Z₂ (g₀² − (2/27) δ g₀²)"], f"Z₂ = 1/{N(V('gemini_Z2_denominator'))}, g₀² = {N(V('gemini_gA_bare_sq'))}, δ = {N(V('gemini_gA_delta'))}",
      [f"From the stated formula = {N(gA_f,4)}; Gemini's claim = {N(V('gemini_gA_claim'))} (needs a hand-added {N(V('gemini_gA_hand_add'))})"],
      f"{N(V('gA_measured'))} — {SRC('gA_measured')}", f"within {N(V('tolerance_gA')*100,2)}%", j_rel(gA_f, V("gA_measured"), V("tolerance_gA")))
mu_c = [Z2 * v for v in V("gemini_mu_bare")]
check("K8c", "C", "Gemini's proton and neutron magnetic moments", (R + ANCH["K8"], "record K8"),
      ["μ = Z₂ × bare value"], f"Z₂ = {N(Z2,4)}; bare values {N(V('gemini_mu_bare')[0])} and {N(V('gemini_mu_bare')[1])}",
      [f"From the stated factor: {N(mu_c[0],4)} and {N(mu_c[1],4)}; Gemini typed {N(V('gemini_mu_typed')[0])} and {N(V('gemini_mu_typed')[1])}, then added a further hand term (claim {N(V('gemini_mu_claim')[0])}, {N(V('gemini_mu_claim')[1])})"],
      "the record states no measured value here", f"the typed value must equal the formula's value within {N(V('tol_percent_1')*100,2)}%",
      j_rel(V("gemini_mu_typed")[0], mu_c[0], V("tol_percent_1")), "Typed numbers, not computed ones.")
check("K8d", "C", "Pressure inside the proton and the D-term", (R + ANCH["K8"], "record K8"),
      ["von Laue: ∫ p(r) r² dr = 0;  D-term"], f"pressure switch at about {N(V('jlab_pressure_switch_fm'))} fm (Jefferson Lab, 2018); D-term bag model {N(V('dterm_bag'))} versus Jefferson Lab estimate {N(V('dterm_jlab'))}",
      [f"Quark-core radius {N(V('core_radius_lo'))} to {N(V('core_radius_hi'))} fm against the switch at {N(V('jlab_pressure_switch_fm'))} fm; D-term ratio {N(V('dterm_jlab')/V('dterm_bag'),3)}"],
      "extractions of the D-term still differ between groups", "not decided", "NOT TESTABLE YET", "Open comparison, not a fail (record K8): the measured value itself is disputed.")
# K9 compactness
Csun = V("compactness_sun"); Cwd = 2 * V("sirius_b_mass") * V("Msun_m") / V("sirius_b_radius_m")
Rns = r14 * 1e3; Cns = 2 * STDM * V("Msun_m") / Rns; Cb = V("buchdahl_limit_num") / V("buchdahl_limit_den")
check("K9", "C", "The containment ladder: compactness C = 2GM/(Rc²)", (R + ANCH["K9"], "record K9"),
      ["C = 2GM/(Rc²);  Buchdahl: C < 8/9 for any static star"],
      f"Sun C = {N(Csun)}; white dwarf (M = {N(V('sirius_b_mass'))} solar masses, R = {N(V('sirius_b_radius_m'))} m); neutron star R = {N(r14,3)} km (stage A3b)",
      [f"White dwarf C = {N(Cwd,3)}; neutron star C = {N(Cns,3)}; Buchdahl limit {fr(V('buchdahl_limit_num'), V('buchdahl_limit_den'))} = {N(Cb,4)}; horizon C = {N(V("horizon_compactness"))}"],
      "Buchdahl (1959): no static star passes the limit without collapsing", "neutron star below the Buchdahl limit and above the white dwarf", "PASS" if Cwd < Cns < Cb else "FAIL")
# K10 Fermi muon lifetime
check("K10b", "C", "Muon lifetime from Fermi theory", (R + ANCH["K10"], "record K10"),
      ["τ_μ = 192π³ℏ / (G_F² m_μ⁵)"], f"G_F = {N(GF)} MeV⁻², m_μ = {N(mmu,7)} MeV",
      [f"τ_μ = {N(tau_mu_c,5)} s", rec_line("K10 muon", tau_mu_c * 1e6, CL("K10mu"), 4)],
      f"{N(V('tau_mu_s'),7)} s — {SRC('tau_mu_s')}", f"within {N(V('tol_percent_1')*100,2)}%", j_rel(tau_mu_c, V("tau_mu_s"), V("tol_percent_1")),
      "The record gives this in microseconds; the comparison above converts it.")

# B8 (new count)
trpg_items = V("inputs_trpg_items"); sm_items = V("inputs_sm_items")
check("B8", "B", "Count of measured inputs", (R + ANCH["P1"], "record P1"),
      ["inputs (TRPG) versus inputs (Standard Model with neutrino masses, gravity and dark energy)"],
      "the record's own lists (equation page and Q2); B and the glueball mass are counted as consequences of the strong force",
      [f"TRPG list holds {I(len(trpg_items))} items; Standard Model list holds {I(len(sm_items))}; saving = {I(len(sm_items) - len(trpg_items))}"],
      f"the record states about {I(CL('inputs_trpg'))} against about {I(CL('inputs_sm'))}", "same counts as the record",
      "PASS" if len(trpg_items) == CL("inputs_trpg") and len(sm_items) == CL("inputs_sm") else "FAIL")
# ---- rev3 B9: B set by the proton (published script, run as published)
pt = PUB["proton_sets_B"][2]
def _first(pat, txt, grp=1, nth=0):
    m_ = re.findall(pat, txt); return float(m_[nth][grp - 1] if isinstance(m_[nth], tuple) else m_[nth])
try:
    B1975 = _first(r"B = ([\d.]+) MeV/fm\^3 \(B\^1/4", pt, nth=0); B_p = _first(r"B = ([\d.]+) MeV/fm\^3 \(B\^1/4", pt, nth=1)
    Rn = _first(r"proton bag radius R = ([\d.]+) fm", pt); rms_p = _first(r"quark-core rms = ([\d.]+) fm", pt)
    cl_lo = _first(r"coupling x0.95: B = ([\d.]+)", pt); cl_hi = _first(r"coupling x1.05: B = ([\d.]+)", pt); parsed_ok = True
except Exception:
    parsed_ok = False; B1975 = B_p = Rn = rms_p = cl_lo = cl_hi = float("nan")
check("B9a", "B", "Containment depth B set by the proton and the Delta, nothing adjusted", (R + ANCH["Q3"], "record Q3"),
      ["fit B and \u03b1_c to M_N and M_\u0394;  E = (4\u03c0/3) B R\u00b3 + ...;  centre-of-mass correction M = \u221a(E\u00b2 \u2212 \u27e8P\u00b2\u27e9);  pion cloud (cloudy bag)"],
      "the published proton_sets_B script, run with python -I (its real output is in stage F2)",
      [f"B = {N(B_p,3)} MeV/fm\u00b3, band {N(cl_lo,3)} to {N(cl_hi,3)} for a \u00b1 coupling change; bag radius {N(Rn,3)} fm; quark-core rms {N(rms_p,3)} fm",
       rec_line("B9 B", B_p, CL("Q3B"), 3), rec_line("B9 band low", cl_lo, CL("Q3Blo"), 3), rec_line("B9 band high", cl_hi, CL("Q3Bhi"), 3), rec_line("B9 rms", rms_p, CL("Q3core"), 3)],
      f"no direct measurement of B; the fit reproduces the measured proton and Delta masses", f"within {N(V('tol_percent_1')*100,2)}% of the record's central value", j_rel(B_p, CL("Q3B"), V("tol_percent_1")) if parsed_ok else "FAIL")
check("B9b", "B", "Check of the same code on the earlier bag fit", (R + ANCH["Q3"], "record Q3"), ["same code with the earlier settings"],
      "settings as printed by the script (no centre-of-mass correction, no cloud)",
      [f"B = {N(B1975,3)} MeV/fm\u00b3", rec_line("B9 earlier fit", B1975, CL("Q3check1975"), 3)],
      f"published about {N(CL('Q3pub'))} MeV/fm\u00b3 (DeGrand, Jaffe, Johnson and Kiskis, the earlier bag fit)",
      f"within {N(V('tol_percent_5')*100,2)}%", j_rel(B1975, CL("Q3pub"), V("tol_percent_5")) if parsed_ok else "FAIL")
# K3d: does the new Q3 radius agree with what K3 says?
check("K3d", "C", "The proton core size: K3's statement against the Q3 run", (R + ANCH["K3"], "record K3"),
      ["compare the quark-core rms of the Q3 run with the range stated in K3"],
      f"Q3 run: rms {N(rms_p,3)} fm; K3 states {N(CL('K3rms_lo'))} to {N(CL('K3rms_hi'))} fm",
      [f"Difference from the upper end of K3's range: {N((rms_p/CL('K3rms_hi')-1)*100,2)}%"], f"{N(CL('K3rms_lo'))} to {N(CL('K3rms_hi'))} fm \u2014 record K3",
      f"within {N(V('tol_percent_5')*100,2)}% of the range", "PASS" if rms_p <= CL("K3rms_hi") * (1 + V("tol_percent_5")) else "PARTIAL",
      "The later Q3 run is the better source; K3's own B figures (stages K3a and K3c) are the older text.")
# Q5: the quark link-pattern idea
qt = PUB["quark_link_pattern"][2]
q5 = re.findall(r"r\^2 = ([\d.]+) \+/- ([\d.]+);\s+simple value ([\d./]+) = ([\d.]+) is (-?[\d.]+) sigma", qt)
q5ok = len(q5) == 2
upr, dnr = (float(q5[0][0]), float(q5[1][0])) if q5ok else (float("nan"), float("nan"))
upsig, dnsig = (float(q5[0][4]), float(q5[1][4])) if q5ok else (float("nan"), float("nan"))
check("Q5", "D", "Quark sectors in the leptons' relational form: the tested idea fails", (R + ANCH["Q5"], "record Q5"),
      ["r\u00b2 = (3Q \u2212 1)/2;  idea tested: r\u00b2 = 3/6, 4/6, 5/6 for leptons, down-type, up-type"],
      "the published quark_link_pattern script, run with python -I (real output in stage F3)",
      [f"up-type r\u00b2 = {N(upr,4)}, which is {N(abs(upsig),3)} standard deviations from the simple value; down-type r\u00b2 = {N(dnr,4)}, {N(abs(dnsig),3)} standard deviations away",
       rec_line("Q5 up", upr, CL("Q5up"), 3), rec_line("Q5 down", dnr, CL("Q5down"), 3)],
      "the measured quark masses (PDG 2024) give these r\u00b2 values", "the simple fractions within 2 standard deviations (look-elsewhere chance in the script's last line)",
      "PASS" if q5ok and abs(upsig) <= 2 and abs(dnsig) <= 2 else "FAIL", "The idea fails and the record drops it. Quark masses stay measured inputs.")

# ================================================================================================ SECTION D: chapter M and P extras
SECTIONS["D"] = "Chapters M and P: the rest"
# M5 wave grid
ng = int(V("wave_n_grid")); xg = np.linspace(0, 2 * np.pi, ng); X, Yg = np.meshgrid(xg, xg)
zx, zy = -np.sin(X), -np.sin(Yg); zxx, zyy = -np.cos(X), -np.cos(Yg)
Kc = (zxx * zyy) / (1 + zx ** 2 + zy ** 2) ** 2; dA = np.sqrt(1 + zx ** 2 + zy ** 2) * (xg[1] - xg[0]) ** 2
gb = float((Kc * dA)[:-1, :-1].sum()); fpos = float((dA * (Kc > 0)).sum() / dA.sum()); fneg = float((dA * (Kc < 0)).sum() / dA.sum())
# critical point classification at the four named points
def classify(x0, y0):
    hxx, hyy = -math.cos(x0), -math.cos(y0); det = hxx * hyy
    return "saddle" if det < 0 else ("peak" if hxx < 0 else "trough")
cps = [classify(0, 0), classify(math.pi, math.pi), classify(math.pi, 0), classify(0, math.pi)]
morse = cps.count("peak") - cps.count("saddle") + cps.count("trough")
check("M5", "D", "Crossing waves: peaks, saddles, troughs and flatness", (R + ANCH["M5"], "record M5"),
      ["z = cos x + cos y;  K = z_xx z_yy / (1 + z_x² + z_y²)²;  ∫ K dA = 0 over one cell (Gauss–Bonnet);  peaks − saddles + troughs = 0 (Morse)"],
      f"grid of {I(ng)}×{I(ng)} points over one cell; critical points at (0,0), (π,π), (π,0), (0,π)",
      [f"Types found: {', '.join(cps)}; peaks − saddles + troughs = {N(morse,3)}", f"Area curved like a sphere = {N(fpos,3)}, like a saddle = {N(fneg,3)}", f"Total curvature over the cell = {N(gb,3)}"],
      f"flat on average: Planck collaboration (2018) finds curvature {N(V('Omega_k'))} ± {N(V('Omega_k_sigma'))}", f"total curvature below {N(V('gauss_bonnet_tol'))} and Morse sum zero",
      "PASS" if abs(gb) < V("gauss_bonnet_tol") and morse == 0 else "FAIL")
# M6 waterfall
def Qsg(vals, signs): rr = [s * math.sqrt(v) for v, s in zip(vals, signs)]; return sum(vals) / sum(rr) ** 2
def solve_light(h2, h1, sgn):
    f = lambda x: Qsg([x, h2, h1], [sgn, 1, 1]) - Q23
    xs = [10 ** (k_ / 200) for k_ in range(-600, 1200)]; out = []
    for a, b in zip(xs, xs[1:]):
        if f(a) * f(b) < 0: out.append(brentq(f, a, b))
    return out
Qcbt = Qsg([qk["c"], qk["b"], qk["t"]], [1, 1, 1]); Qscb = Qsg([qk["s"], qk["c"], qk["b"]], [-1, 1, 1])
chs = solve_light(qk["b"], qk["t"], 1); cc_ = min(chs, key=lambda v: abs(v - qk["c"])); sst = solve_light(cc_, qk["b"], -1)
ss_ = min(sst, key=lambda v: abs(v - qk["s"])) if sst else float("nan")
wst = V("waterfall_stated")
check("M6a", "D", "Quark masses from linked Koide triples: charm", (R + ANCH["M6"], "record M6"),
      ["Q(c, b, t) = 2/3 solved for the charm mass, starting from b and t only"], f"m_b = {N(qk['b'])} MeV, m_t = {N(qk['t'])} MeV",
      [f"Q(c,b,t) = {N(Qcbt,4)}; predicted charm = {N(cc_,5)} MeV; measured charm = {N(qk['c'])} MeV ({N((cc_/qk['c']-1)*100,2)}% off)", rec_line("M6 charm", cc_, wst["charm_pred"], 4)],
      f"{N(qk['c'])} MeV — {SRC('q_c')}", f"within {N(V('tol_percent_1')*100,2)}%", j_rel(cc_, qk["c"], V("tol_percent_1")),
      "Quark masses run with energy and these inputs mix scales, so the record calls this a hint.")
check("M6b", "D", "Quark masses from linked Koide triples: strange", (R + ANCH["M6"], "record M6"),
      ["Q(−s, c, b) = 2/3 solved for the strange mass"], f"m_c (from the step above) = {N(cc_,5)} MeV, m_b = {N(qk['b'])} MeV",
      [f"Q(s,c,b) with the measured values = {N(Qscb,4)}; predicted strange = {N(ss_,4)} MeV; measured {N(qk['s'])} MeV ({N((ss_/qk['s']-1)*100,2)}% off)", rec_line("M6 strange", ss_, wst["strange_pred"], 4)],
      f"{N(qk['s'])} MeV — {SRC('q_s')}", f"within {N(V('tol_percent_1')*100,2)}%", j_rel(ss_, qk["s"], V("tol_percent_1")))
# M7 resistance
rM = V("Mars_orbit_in_AU") * V("AU_m"); rE = V("AU_m")
shap = (2 * G * MSUN / C ** 3) * math.log(4 * rE * rM / RSUN ** 2) * 1e6
rt = 2 * shap; lo_rt, hi_rt = V("shapiro_roundtrip_range_us")
check("M7a", "D", "Resistance as refractive index: radar delay grazing the Sun", (R + ANCH["M7"], "record M7"),
      ["Δt = (2GM/c³) ln(4 r_E r_M / b²)  (one way)"], f"r_E = {N(rE)} m, r_M = {N(V('Mars_orbit_in_AU'))} AU, b = R_sun = {N(RSUN)} m",
      [f"One way = {N(shap,3)} microseconds; round trip = {N(rt,3)} microseconds", rec_line("M7a one way", shap, CL("M7a"), 3)],
      f"{N(lo_rt)} to {N(hi_rt)} microseconds round trip — {SRC('shapiro_roundtrip_range_us')}", "round trip inside the quoted range", j_range(rt, lo_rt, hi_rt))
# M7 bands & localization
def tr(k_, n1=1.0, n2=1.5, a=0.5, bb=0.5):
    def layer(n, d): q = k_ * n; return np.array([[math.cos(q * d), math.sin(q * d) / n], [-n * math.sin(q * d), math.cos(q * d)]])
    T = layer(n1, a) @ layer(n2, bb); return abs(np.trace(T)) / 2
ks = np.linspace(0.01, 12, 4000); allowed = np.array([tr(k_) <= 1 for k_ in ks]); edges = np.where(np.diff(allowed.astype(int)) != 0)[0]
gaps = [(ks[edges[i]], ks[edges[i + 1]]) for i in range(0, len(edges) - 1, 2) if not allowed[edges[i] + 1]]
check("M7b", "D", "Periodic resistance makes bands and gaps", (R + ANCH["M7"], "record M7"),
      ["|Tr T|/2 ≤ 1 allowed;  > 1 forbidden (transfer matrix of alternating layers)"], "alternating layers of two refractive indices, wave number scanned from small to large",
      [f"Forbidden gaps found: {', '.join(f'{N(a_,3)} to {N(b_,3)}' for a_, b_ in gaps[:6])}"],
      "Bloch (1928); Kronig and Penney (1931): periodic media have forbidden gaps", f"at least {I(V('bands_min_gaps'))} gaps", "PASS" if len(gaps) >= V("bands_min_gaps") else "FAIL")
rng3 = np.random.default_rng(0)
def T_random(Nl, k_=3.0, spread=0.3):
    T = np.eye(2)
    for _ in range(Nl):
        n = 1 + spread * rng3.random(); d = 0.5; q = k_ * n
        T = np.array([[math.cos(q * d), math.sin(q * d) / n], [-n * math.sin(q * d), math.cos(q * d)]]) @ T
    t = 2 / abs(T[0, 0] + T[1, 1] + 1j * (T[0, 1] * 1 - T[1, 0] / 1)); return t * t
trans = [(Nl, float(np.mean([T_random(Nl) for _ in range(20)]))) for Nl in V("resistance_layers")]
check("M7c", "D", "Random resistance traps waves (Anderson localization)", (R + ANCH["M7"], "record M7"),
      ["transmitted fraction through N random layers falls exponentially with N"], "layer counts " + ", ".join(I(n) for n, _ in trans),
      ["Transmitted fraction: " + "; ".join(f"{I(n)} layers → {N(t_,3)}" for n, t_ in trans)],
      "Anderson (1958)", f"fall by at least a factor {I(V('localization_factor'))} from the shortest to the longest stack", "PASS" if trans[0][1] / trans[-1][1] >= V("localization_factor") else "FAIL")
# M8 compactness of Hubble sphere
Ch = 2 * G * Mu / (Rh * V("c_exact") ** 2)
check("M8", "D", "The universe on the containment ladder", (R + ANCH["M8"], "record M8"),
      ["C = 2GM/(Rc²) with R = c/H₀ and M = (4π/3)R³ρ_c"], f"R = {N(Rh,3)} m, M = {N(Mu,3)} kg",
      [f"C = {N(Ch,5)}"], "identity for a flat universe", f"equal to 1 within {N(V('tol_lepton_mass')*100,2)}%", j_rel(Ch, 1.0, V("tol_lepton_mass")),
      "An identity of flat cosmology: consistent with the idea, not a proof of it.")
# M9 factor 2
full = 4 * G * MSUN / (RSUN * C * C) * ASEC
cases = [("scalar resistance (Nordström 1913)", None), ("time-only resistance (Newton, Soldner 1801)", 0.0), ("time and space in equal shares (Einstein 1915)", 1.0)]
for n_, (nm, gam) in enumerate(cases):
    dv = 0.0 if gam is None else (1 + gam) / 2 * full
    check(f"M9-{'abc'[n_]}", "D", f"Light bending with {nm}", (R + ANCH["M9"], "record M9"),
          ["δ = (1 + γ)/2 × 4GM/(Rc²)"], ("scalar: no coupling to light, so zero" if gam is None else f"γ = {N(gam)}") + f"; 4GM/(Rc²) = {N(full,4)} arcsec",
          [f"Deflection = {N(dv,4)} arcsec"], f"{N(V('bending_measured'))} arcsec — {SRC('bending_measured')}", f"within {N(V('tol_percent_1')*100,2)}%", j_rel(dv, V("bending_measured"), V("tol_percent_1")),
          "Ruled out by the data." if j_rel(dv, V("bending_measured"), V("tol_percent_1")) == "FAIL" else f"The only case that survives. Cassini, as the record cites it: \u03b3 \u2212 1 = {N(V('gamma_minus_1'))} \u00b1 {N(V('gamma_minus_1_sigma'))}.")
# M10 reach
reachP = V("hc_eVm") / V("p_ripple_eV"); reachG = V("hc_eVm") / V("graviton_limit_eV"); need = V("gravity_reach_needed_m")
check("M10a", "D", "The containment field cannot be gravity: its reach", (R + ANCH["M10"], "record M10"),
      ["reach ≈ hc / (mc²)"], f"P ripple mass = {N(V('p_ripple_eV'))} eV (the glueball); graviton limit = {N(V('graviton_limit_eV'))} eV (LIGO-Virgo-KAGRA 2021)",
      [f"Reach of P = {N(reachP,3)} m ({N(reachP*1e15,3)} fm); reach allowed for gravity = at least {N(reachG,3)} m"],
      f"gravity acts over cosmic distances: {SRC('graviton_limit_eV')}", f"for P to be gravity its reach would have to be at least {N(need)} m", "PASS" if reachP < need else "FAIL",
      "PASS here means the check correctly shows P is not gravity.")
check("M10b", "D", "All layers must feel gravity alike (MICROSCOPE)", (R + ANCH["M10"], "record M10"),
      ["Δa/a < measured bound for every pair of materials"], "TRPG layer-dependent gravity has not been computed",
      ["No TRPG number exists yet to compare with the measured bound."], f"MICROSCOPE team (2022) as the record cites", "a TRPG value of the layer difference", "NOT TESTABLE YET")
# M11 expansion
sh = lambda a: ((V("Omega_m") / a ** 3) / (Hf(a) / V("H0")) ** 2, (V("Omega_r") / a ** 4) / (Hf(a) / V("H0")) ** 2, (OLs) / (Hf(a) / V("H0")) ** 2)
a_list = V("expansion_sizes")
lines = []
for a in a_list:
    m_, r2, l_ = sh(a); lines.append(f"size {N(a,1)} of today: matter {N(m_*100,3)}%, radiation {N(r2*100,3)}%, dark energy {N(l_*100,3)}%")
tot_ok = all(abs(sum(sh(a)) - 1) < V("tol_sum_one") for a in a_list)
check("M11", "D", "Always full: the contents add to the whole at every size", (R + ANCH["M11"], "record M11"),
      ["Ω_m + Ω_r + Ω_Λ = 1;  shares at size a: Ω_i(a) / E(a)²"], f"Ω_m = {N(V('Omega_m'))}, Ω_r = {N(V('Omega_r'))}, Ω_Λ = {N(OLs,5)}",
      lines + [f"Shares sum to {I(V('percent'))}% at each size: {'yes' if tot_ok else 'no'}"],
      f"Planck 2018: curvature {N(V('Omega_k'))} ± {N(V('Omega_k_sigma'))}", f"sum equal to 1 within {N(V('tol_sum_one'))} and curvature consistent with zero", "PASS" if tot_ok and abs(V("Omega_k")) <= V("Omega_k_sigma") * 1.0 + 1e-12 + 0.5 * V("Omega_k_sigma") else "FAIL",
      "The curvature term is the only 'empty' share left over; Planck's value is consistent with zero.")
# M12 uncertainty
gp = V("gauss_pulse"); x = np.linspace(-gp["extent"], gp["extent"], gp["points"]); s0 = gp["width"]; psi = np.exp(-x ** 2 / (4 * s0 * s0)); pr = abs(psi) ** 2; pr /= np.trapezoid(pr, x)
dx = math.sqrt(np.trapezoid(x ** 2 * pr, x)); kk = np.fft.fftshift(np.fft.fftfreq(len(x), x[1] - x[0])) * 2 * np.pi
Pk = abs(np.fft.fftshift(np.fft.fft(psi))) ** 2; Pk /= np.trapezoid(Pk, kk); dk = math.sqrt(np.trapezoid(kk ** 2 * Pk, kk))
check("M12", "D", "Uncertainty is wave math", (R + ANCH["M12"], "record M12"),
      ["Δx Δk ≥ 1/2  (equality for a Gaussian pulse)"], f"Gaussian pulse of width {N(s0)} on {N(gp['points'])} points",
      [f"Δx Δk = {N(dx*dk,5)}"], "Heisenberg (1927): the minimum is one half", f"within {N(V('tol_lepton_mass')*100,2)}% of one half", j_rel(dx * dk, 0.5, V("tol_lepton_mass")))
# M14 entanglement density
lp = V("planck_length_m"); dens = 1 / (4 * lp ** 2)
check("M14", "D", "Entanglement per square metre of space", (R + ANCH["M14"], "record M14"),
      ["one unit per 4 Planck areas:  1 / (4 ℓ_P²)"], f"ℓ_P = {N(lp)} m",
      [f"Density = {N(dens,3)} per m²", rec_line("M14 density", dens, V("entangle_density_stated"), 3)],
      f"the record states about {N(V('entangle_density_stated'))} per square metre", f"within a factor {I(V('factor_three_tol'))} of the stated figure",
      "PASS" if (1 / V("factor_three_tol")) <= dens / V("entangle_density_stated") <= V("factor_three_tol") else "FAIL",
      f"The record's figure of {N(V('entangle_density_stated'))} per square metre is one power of ten below what the arithmetic gives. The weak-gravity argument does not depend on it; the figure in the record is a slip.")
# M15 small world
sw = V("small_world_cases"); swS = V("small_world_stated")
lines = []; okw = True
for (Nn, kn), st in zip(sw, swS):
    v_ = math.log(Nn) / math.log(kn); lines.append(f"N = {N(Nn,2)}, {N(kn,3)} links each → {N(v_,3)} steps (record: {N(st,3)})"); okw &= abs(v_ / st - 1) <= V("tol_percent_5")
check("M15a", "D", "Small-world estimate of steps between any two things", (R + ANCH["M15"], "record M15"),
      ["steps ≈ ln N / ln k"], f"people ({N(sw[0][0],2)} with {N(sw[0][1],3)} close contacts) and the universe's particles ({N(sw[1][0],2)} with {N(sw[1][1],3)} links each, and with {N(sw[2][1],3)} links each)", lines,
      "Watts and Strogatz (1998)", f"within {N(V('tol_percent_5')*100,2)}% of the record's figures", "PASS" if okw else "FAIL")
# M16 area law numbers
check("M16", "D", "Area law and the strength of gravity", (R + ANCH["M16"], "record M16"),
      ["same test as stage B1a"], "see stage B1a", [f"Boundary-fit error {N(rb,3)} versus volume-fit error {N(rv,3)}"],
      "Bombelli, Koul, Lee and Sorkin (1986); Srednicki (1993)", f"ratio at least {N(V('area_ratio_min'))}", "PASS" if rv / rb >= V("area_ratio_min") else "FAIL")
# M18 balance angles
def Qang(v): Qv = sum(v) / sum(math.sqrt(t) for t in v) ** 2; return Qv, math.degrees(math.acos(math.sqrt(1 / (3 * Qv))))
a_l = Qang([me, mmu, mtau])[1]; a_d = Qang([qk["d"], qk["s"], qk["b"]])[1]; a_u = Qang([qk["u"], qk["c"], qk["t"]])[1]
qas = V("quark_angle_stated")
check("M18a", "D", "The balance angle: leptons", (R + ANCH["M18"], "record M18"),
      ["θ = arccos √(1/(3Q))  (45° when Q = 2/3)"], "charged-lepton masses as in stage A4", [f"Angle = {N(a_l,4)} degrees", rec_line("M18 lepton angle", a_l, qas["lepton"], 4)],
      "45° — Koide (1981); Foot (1994)", f"within {N(V('tol_percent_1')*100,2)}%", j_rel(a_l, math.degrees(math.acos(math.sqrt(1 / (3 * Q23)))), V("tol_percent_1")))
a_ref = math.degrees(math.acos(math.sqrt(1 / (3 * Q23))))
for nm, av_ in (("down-type quarks", a_d), ("up-type quarks", a_u)):
    check("M18-" + nm[:2], "D", f"The balance angle: {nm}", (R + ANCH["M18"], "record M18"),
          ["same formula on the quoted masses"], "quark masses as listed in the inputs", [f"Angle = {N(av_,4)} degrees; lepton angle = {N(a_ref,4)} degrees", rec_line("M18 quark angle " + nm, av_, qas["down" if nm.startswith("down") else "up"], 4)],
          f"the lepton value {N(a_ref,4)} degrees", f"within {N(V('tol_percent_1')*100,2)}% of the lepton angle", j_rel(av_, a_ref, V("tol_percent_1")),
          "Quarks do not sit at the lepton balance. The record calls this unsolved, as in all of physics.")
# P2 quark running
MZ = V("mu_quark_run")["MZ"]; aZ = V("mu_quark_run")["alpha_s_MZ"] / math.pi; mc_, mb_, mt_ = V("mu_quark_run")["mc"], V("mu_quark_run")["mb"], V("mu_quark_run")["mt"]
def b_(nf): return (11 - 2 * nf / 3) / 4, (102 - 38 * nf / 3) / 16
def g_(nf): return 1.0, (202 / 3 - 20 * nf / 9) / 16
def nf_at(mu): return 6 if mu > mt_ else 5 if mu > mb_ else 4 if mu > mc_ else 3
def run(a0, mu0, mu1, lnm=0.0):
    th = [mc_, mb_, mt_]; pts = sorted([mu0, mu1] + [t for t in th if min(mu0, mu1) < t < max(mu0, mu1)], reverse=(mu1 < mu0)); a, l = a0, lnm
    for x0, y0 in zip(pts, pts[1:]):
        nf = nf_at(math.sqrt(x0 * y0)); b0, b1 = b_(nf); g0, g1 = g_(nf)
        f = lambda t, s: [-b0 * s[0] ** 2 - b1 * s[0] ** 3, -g0 * s[0] - g1 * s[0] ** 2]
        a, l = solve_ivp(f, [math.log(x0 * x0), math.log(y0 * y0)], [a, l], rtol=1e-10, atol=1e-12).y[:, -1]
    return a, l
a_at = lambda mu: run(aZ, MZ, mu)[0]
m_to = lambda m0, mu0, mu1: m0 * math.exp(run(a_at(mu0), mu0, mu1)[1])
qi = V("quark_run_inputs"); lines = []; angs = []
for mu in V("quark_run_scales"):
    m = {k_: m_to(m0, mu0, mu) for k_, (m0, mu0) in qi.items()}
    up = Qang([m["u"], m["c"], m["t"]]); dn = Qang([m["d"], m["s"], m["b"]]); angs.append((mu, up[1], dn[1]))
    lines.append(f"at {N(mu,5)} GeV: up-type angle {N(up[1],4)}°, down-type angle {N(dn[1],4)}°")
spread_up = max(a[1] for a in angs) - min(a[1] for a in angs); spread_dn = max(a[2] for a in angs) - min(a[2] for a in angs)
rs = V("quark_angle_run_stated"); at91 = [a for a in angs if abs(a[0] - MZ) < 1e-6][0]
check("P2a", "D", "Quark balance angles are the same at every energy", (R + ANCH["P2"], "record P2"),
      ["two-loop running: dα/d lnμ² = −b₀α² − b₁α³;  d ln m/d lnμ² = −γ₀α − γ₁α²"],
      f"α_s(M_Z) = {N(V('mu_quark_run')['alpha_s_MZ'])}; quark masses at their quoted scales ({SRC('quark_run_inputs')}); energies {', '.join(N(v,5) for v in V('quark_run_scales'))} GeV",
      lines + [f"Spread of the up-type angle across energies = {N(spread_up,3)}°; down-type = {N(spread_dn,3)}°", rec_line("P2 up angle at M_Z", at91[1], rs["up"], 4), rec_line("P2 down angle at M_Z", at91[2], rs["down"], 4)],
      "the angle should not change with energy because QCD rescales all quark masses by one factor (record P2)", f"spread below {N(V('tol_percent_1'))} degrees",
      "PASS" if spread_up < V("tol_percent_1") and spread_dn < V("tol_percent_1") else "FAIL")
check("P2b", "D", "Do the quarks sit at the lepton balance once run to one energy?", (R + ANCH["P2"], "record P2"),
      ["compare the angles with the lepton angle"], f"angles at {N(MZ)} GeV: up-type {N(at91[1],4)}°, down-type {N(at91[2],4)}°; leptons {N(a_ref,4)}°",
      [f"Up-type is {N(at91[1]-a_ref,3)}° and down-type {N(at91[2]-a_ref,3)}° away from the leptons."], f"the lepton value {N(a_ref,4)} degrees",
      f"within {N(V('tol_percent_1')*100,2)}% of the lepton angle", "PASS" if abs(at91[1] / a_ref - 1) <= V("tol_percent_1") and abs(at91[2] / a_ref - 1) <= V("tol_percent_1") else "FAIL",
      "They do not. The record states this too: it is now a sharp target for TRPG's curvature-phase idea, not a result.")
# P4 gravity tightened
Lt = np.array(V("area_blocks_tidy"))
S_cross = ent(Ng, mg, True, Lt, V("crossing_amplitude"))
eb = float(np.sqrt(np.mean((np.polyval(np.polyfit(4 * Lt, S_cross, 1), 4 * Lt) - S_cross) ** 2))); ev_ = float(np.sqrt(np.mean((np.polyval(np.polyfit(Lt * Lt, S_cross, 1), Lt * Lt) - S_cross) ** 2)))
check("P4a", "D", "The area law in the crossing-wave medium", (R + ANCH["P4"], "record P4"),
      ["couplings k = 1 + a cos x cos y;  fit S to boundary and to volume"], f"a = {N(V('crossing_amplitude'))}; blocks of side {', '.join(I(l) for l in Lt)}",
      [f"Boundary-fit error {N(eb,3)}; volume-fit error {N(ev_,3)}; ratio {N(ev_/eb,3)}", rec_line("P4 ratio", ev_ / eb, CL("P4"), 3)],
      "Bombelli, Koul, Lee and Sorkin (1986); Srednicki (1993)", f"boundary fit better by at least a factor {I(V('area_ratio_min_crossing'))}", "PASS" if ev_ / eb >= V("area_ratio_min_crossing") else "FAIL")
coefs = []
for m in V("layer_masses"):
    Sm_ = ent(Ng, m, False, Lt); coefs.append((m, float(np.polyfit(4 * Lt, Sm_, 1)[0])))
check("P4b", "D", "Every layer contributes to G, and G acts on all of them alike", (R + ANCH["P4"], "record P4"),
      ["area coefficient c(m) for fields of different mass; G is the single total of the contributions (Jacobson 1994; Susskind and Uglum 1994)"],
      "layer masses " + ", ".join(N(m) for m, _ in coefs), ["Coefficient per boundary unit: " + "; ".join(f"mass {N(m)} → {N(c_,3)}" for m, c_ in coefs)],
      "no measured value: this stage only shows layers differ in what they give to G", "coefficients decrease with mass", "PASS" if all(coefs[i][1] > coefs[i + 1][1] for i in range(len(coefs) - 1)) else "FAIL",
      "It does not test the equivalence principle itself (MICROSCOPE bound, stage M10b).")
# P5 twist
def delta_of(a, b, c_): rr = np.sqrt([a, b, c_]); M_ = rr.mean(); return math.acos((rr[2] / M_ - 1) / math.sqrt(2))
d0 = delta_of(me, mmu, mtau); rng4 = np.random.default_rng(int(V("twist_seed")))
ds = [delta_of(me, mmu + rng4.normal(0, V("sigma_m_mu")), mtau + rng4.normal(0, V("sigma_m_tau"))) for _ in range(int(V("twist_samples")))]
sig = float(np.std(ds)); nsig = (d0 - D29) / sig; ts = V("twist_stated")
check("P5", "D", "The twist δ = 2/9 from the data", (R + ANCH["P5"], "record P5"),
      ["δ = arccos((√m_τ / M − 1)/√2),  M = mean of the three square roots;  compare 2/9"],
      f"masses and uncertainties as in the inputs; {I(V('twist_samples'))} random draws",
      [f"δ from data = {N(d0,6)} ± {N(sig,2)} rad; 2/9 = {N(D29,6)}; difference = {N(nsig,2)} standard deviations", rec_line("P5 twist", d0, ts[0], 6)],
      f"2/9 — Brannen (2006)", f"within {I(V('twist_sigma_limit'))} standard deviations", "PASS" if abs(nsig) <= V("twist_sigma_limit") else "FAIL",
      "That 2/9 equals (2/3)/3 is the project's conjecture, not derived (record P5).")

# ================================================================================================ SECTION E: not recomputed here
SECTIONS["E"] = "Stated in the record but not recomputed on this page"
NOTRERUN = [
 ("J1/G1", "The band test on the massive elementary particles (record J1, G1: a weak hint, p about a tenth)", "no script for this test is among the downloads yet"),
 ("G1", "The golden-ratio ladder test", "no script among the downloads yet"),
 ("M4/M13/M17", "The scorecards (the record's own PASS and PARTIAL counts)", "these are the record's own judgements; the stages above are the evidence they rest on"),
 ("P1", "The comparison count for the Standard Model plus gravity and dark energy", "not recounted here"),
 ("K8", "That the neutron-minus-proton mass and the von Laue residual are built in by construction", "this is a reading of Gemini's code, not a number"),
]
for cid, ttl, why in NOTRERUN:
    check("X-" + cid.replace("/", "-"), "E", ttl, (R + "part-1-every-idea-step-by-step", "record"), ["—"], "not run", [why], "as the record states", "not applicable", "NOT RE-RUN")

# ================================================================================================ SECTION F: the published scripts, run as published
SECTIONS["F"] = "The published scripts, run as published"
for n_, pre_ in enumerate(PUBLISHED):
    name, rc, out_, err_ = PUB[pre_]
    miss = re.findall(r"No module named '(\w+)'", err_)
    ran = rc == 0 and out_.strip() != ""
    tail = err_.strip().splitlines()[-1] if err_.strip() else ""
    check(f"F{n_ + 1}", "F", f"{name}", (R + "chapter-q-the-line-by-line-check-and-the-neutron-star-prediction-14-nov-evening", "record Chapter Q"),
          [f"python -I {name}   (working directory separate from the downloads)"],
          "run with the isolated interpreter in a separate working directory, a copy of the workbook beside it",
          ([f"Exit status {I(rc)}: ran to the end."] if ran else [f"Exit status {I(rc if rc is not None else 0)}: DID NOT RUN HERE. {tail}"]),
          "the script's own printed output follows, unedited", "the script must run to the end in this environment",
          "PASS" if ran else "FAIL",
          ("" if ran else (f"The module {', '.join(miss)} is not installed in this environment, so this script is not verified here. The page computes the same result itself in stages B5 and K10a." if miss else "")),
          pre=(out_.strip() + (("\n--- error output ---\n" + err_.strip()[-1200:]) if err_.strip() else "")))

# ================================================================================================ render
def esc(s): return html.escape(s, quote=False)
def badge(st): return f'<span class="badge badge-{CLS[st]}">{st}</span>'
counts = {s: sum(1 for c in CHECKS if c["status"] == s) for s in STATUSES}
def render_check(c):
    out = [f'<article class="check" id="c-{esc(c["id"])}">',
           f'<header><h3>{esc(c["id"])}. {esc(c["title"])}</h3>{badge(c["status"])}</header>']
    out += [f'<div class="eq">{esc(e)}</div>' for e in c["eqs"]]
    out.append('<dl class="calc">')
    out.append(f'<dt>Numbers in</dt><dd>{esc(c["subst"])}</dd>')
    out.append('<dt>Computed</dt><dd>' + "<br>".join(esc(l) for l in c["computed"]) + '</dd>')
    out.append(f'<dt>Measured</dt><dd>{esc(c["measured"])}</dd>')
    out.append(f'<dt>Tolerance</dt><dd>{esc(c["tol"])}</dd>')
    out.append(f'<dt>Result</dt><dd>{badge(c["status"])}' + (f' {esc(c["note"])}' if c["note"] else "") + '</dd>')
    out.append('</dl>')
    if c.get('pre'): out.append('<pre class="src"><code>' + esc(c['pre']) + '</code></pre>')
    out.append(f'<p class="note">In the record: <a href="{c["rec"][0]}">{esc(c["rec"][1])}</a></p></article>')
    return "\n".join(out)

toc = ["<ul>"]
for key, name in SECTIONS.items():
    toc.append(f'<li><a href="#sec-{key.lower()}">{esc(name)}</a>: ' + ", ".join(f'<a href="#c-{esc(c["id"])}">{esc(c["id"])}</a> {badge(c["status"])}' for c in CHECKS if c["sec"] == key) + "</li>")
toc.append("</ul>")
summary_rows = "".join(f'<tr><td>{badge(s)}</td><td>{I(counts[s])}</td></tr>' for s in STATUSES)
body = [
 "<h1>The math, step by step</h1>",
 '<p class="lede">Every equation on this site, written out, with the numbers put in, the answer, the measured value, and a verdict you can check.</p>',
 f'<p>This page was written by a program, not by hand. The program reads the measured values (each with its source) from one data file, does every calculation itself, and then decides PASS or FAIL by comparing with a tolerance that is printed next to the answer. No number on this page is typed in. Where the record says a result failed or is only partly there, this page shows FAIL or PARTIAL. The same code is on the <a href="code.html">code page</a>, so you can copy it and run it yourself.</p>',
 f'<p>How to read a verdict. <strong>PASS</strong>: the computed value is inside the stated tolerance. <strong>PARTIAL</strong>: outside it but within {N(PF,2)} times the tolerance, or a hint the record itself calls partial. <strong>FAIL</strong>: outside it. <strong>NOT TESTABLE YET</strong>: no measurement or no TRPG number exists to compare. <strong>INPUT</strong>: a measured value used as an input, not a prediction. <strong>INHERITED</strong>: taken unchanged from established physics. <strong>NOT RE-RUN</strong>: stated in the record but not recomputed here.</p>',
 f'<p>The tally is {I(len(CHECKS))} stages in all: ' + ", ".join(f"{I(counts[s])} {s}" for s in STATUSES if counts[s]) + f'. A verdict applies to one stage only. A PASS on a stage does not make a theory true, and PASS counts here are not the record\'s scorecard (see <a href="record.html#{ANCH["M17"]}">record M17</a>).</p>',
 f'<p>Where the record states a number and this run gives a different one, the stage says so. Number of stages where that happened: {I(len(MISMATCH))}.</p>',
 '<table class="summary"><thead><tr><th>Verdict</th><th>Stages</th></tr></thead><tbody>' + summary_rows + '</tbody></table>',
 '<nav class="toc" aria-label="Stages"><h2>Jump to a stage</h2>' + "\n".join(toc) + '</nav>',
]
for key, name in SECTIONS.items():
    body.append(f'<h2 id="sec-{key.lower()}">{esc(name)}</h2>')
    body += [render_check(c) for c in CHECKS if c["sec"] == key]
body_html = "\n".join(body)
available = {p.name for p in BUILD.glob("*.html")} | {"math.html"}
out = sc.page("math.html", "The math, step by step",
              "Every equation on the TRPG site with the numbers substituted, the computed result, the measured value, the tolerance and a PASS or FAIL for each stage, all computed by the published code.",
              body_html, available)

# ------------------------------------------------------------------------------------------------ provenance check
def main_text(page_html):
    m = re.search(r"<main[^>]*>(.*)</main>", page_html, re.S).group(1)
    m = re.sub(r"<[^>]+>", " ", m)
    return html.unescape(m)
def provenance(page_html):
    bad = []
    for tok in NUM.findall(main_text(page_html)):
        if tok not in EMITTED and tok not in FORMULA: bad.append(tok)
    return bad
sb = static_check()
if sb:
    print("STATIC CHECK FAILED:", sb[:10]); sys.exit(4)
bad = provenance(out)
if bad:
    mt = main_text(out)
    print("PROVENANCE CHECK FAILED: numbers not emitted by the generator:", sorted(set(bad)))
    for tok in sorted(set(bad)):
        mm_ = re.search(r"(?<![A-Za-z_\d.])" + re.escape(tok) + r"(?!\d|\.\d)", mt); i = mm_.start() if mm_ else 0; print("  ", tok, "->", ascii(mt[max(0, i - 50):i + 30]))
    sys.exit(2)
if "--plant" in sys.argv:
    planted = out.replace("</main>", "<p>A typed number " + "".join(["7"] * 4) + " slipped in.</p></main>")
    if not provenance(planted):
        print("PLANT NOT CAUGHT: provenance check is blind"); sys.exit(3)
    print("plant caught: the provenance check rejected a typed number as it must:", sorted(set(provenance(planted))))
(BUILD / "math.html").write_text(out, encoding="utf-8")
sc.rewrite_navs(BUILD); sc.update_sitemap(BUILD, ["math.html"])
print(f"math.html written: {len(CHECKS)} stages; " + ", ".join(f"{counts[s]} {s}" for s in STATUSES if counts[s]))
print("record mismatches:", MISMATCH)
print(f"provenance: {len(set(NUM.findall(main_text(out))))} distinct numbers on the page, all emitted by the generator or part of a printed formula")