#!/usr/bin/env python3
"""
verify_P059.py -- Addendum 59: OP-B-Mass top/charm mass-ratio bound.

This verifier checks the numerical situation and the route-closing arithmetic
in 59_Addendum_ConjBMassRatio.tex.  The main conservative conclusion is sound:
Z3 charge alone does not determine mt/mc, and OP-B-Mass is bounded but open.
The flagged issues are stale arithmetic in the route tests: Tr(Xsec^2) and
16*pi^6 are printed about ten times too large, the cube-root trace is off, and
the real moment index where mu_k=1 is about 528, not about 189.
"""

from __future__ import annotations

import math
import sys
from pathlib import Path



PASS = FAIL = 0
_N = 0

def record(label, ok, computed="", claimed="", detail=""):
    """Modern-format check line; behavior-preserving port of verify_common."""
    global PASS, FAIL, _N
    _N += 1
    ok = bool(ok)
    desc = label
    if ok:
        PASS += 1
    else:
        FAIL += 1
        if "Expected" in detail:
            i = detail.find("Expected")
            desc = f"{label} -- {detail[i:]}"
            detail = detail[:i].rstrip().rstrip(";")
    print(f"  [{'PASS' if ok else 'FAIL'}] {_N:>2}. {desc}")
    if computed != "" or claimed != "":
        print(f"        computed: {computed}")
        print(f"        claimed : {claimed}")
    if detail:
        print(f"        {detail}")
    return ok

def check(label, computed, claimed, *, rel=1e-3, abs_tol=None, detail=""):
    if abs_tol is not None:
        ok = abs(computed - claimed) <= abs_tol
        err_detail = f"abs err={abs(computed - claimed):.6g}, tol={abs_tol:.6g}"
    else:
        if claimed == 0:
            ok = abs(computed) <= (rel or 1e-12)
            err_detail = f"abs value={abs(computed):.6g}, tol={rel:.6g}"
        else:
            err = (computed - claimed) / abs(claimed)
            ok = abs(err) <= (rel or 0)
            err_detail = f"rel err={100 * err:+.6g}%, tol={100 * (rel or 0):.6g}%"
    return record(label, ok, computed, claimed, err_detail + (f"; {detail}" if detail else ""))

print("P059 -- OP-B-Mass Top/Charm Bound")

ROOT = Path(__file__).resolve().parents[1]
TEX = (ROOT / "59_Addendum_ConjBMassRatio.tex").read_text()

PI = math.pi
MU0 = 4 * PI**3 + PI**2 + PI
MU1 = 16 * PI**3 / 5 + 3 * PI**2 / 4 + 2 * PI / 3
MU = MU1 / MU0
E_EDGE = PI
E_BOUNDARY = PI**2
E_BULK = 4 * PI**3
E_CHARM = E_EDGE + E_BOUNDARY
GAP = math.log(MU0) / MU
E_TOP = E_CHARM + GAP
ME_GEV = 0.000511
MT = 172.69
MC_2GEV = 1.275
MC_MC = 1.260
MC_POLE = 1.67


def pct(value: float, target: float) -> float:
    return 100 * (value - target) / abs(target)


def mu_k(k: float) -> float:
    return 16 * PI**3 / (k + 4) + 3 * PI**2 / (k + 3) + 2 * PI / (k + 2)


def root_mu_equals_one() -> float:
    lo = 0.0
    hi = 2000.0
    for _ in range(200):
        mid = (lo + hi) / 2
        if mu_k(mid) > 1:
            lo = mid
        else:
            hi = mid
    return (lo + hi) / 2


check("mu0", MU0, 137.036304, rel=5e-9)
check("mu1", MU1, 108.716684, rel=5e-9)
check("MU", MU, 0.793342, rel=5e-7)
check("charm energy Ec", E_CHARM, 13.011, rel=2e-5)
check("gap ln(mu0)/MU", GAP, 6.202, rel=2e-5)
check("predicted top energy", E_TOP, 19.213, rel=7e-6)
check(
    "printed mass-formula exponent Et-Ee",
    E_TOP - E_EDGE,
    16.101,
    rel=5e-4,
    detail="Expected fail: exact Et-Ee is about 16.0715, not 16.101.",
)
check("predicted top mass from mass formula", ME_GEV * math.exp(MU * (E_TOP - E_EDGE)), 176.1, rel=3e-4)

ratio_2gev = MT / MC_2GEV
ratio_mc = MT / MC_MC
ratio_pole = MT / MC_POLE
check("mt/mc at 2 GeV", ratio_2gev, 135.44, rel=3e-5)
check("mu0 discrepancy vs 2 GeV ratio", pct(MU0, ratio_2gev), 1.17, rel=6e-3)
check("mt/mc at mc scale", ratio_mc, 137.1, rel=4e-4)
check("mc-scale residual percent", abs(pct(ratio_mc, MU0)), 0.03, rel=6e-1)
check("mt/mc pole row", ratio_pole, 103.4, rel=8e-5)

check("Z3 top/charm charge difference", (2 - 0) % 3, 2, rel=1e-12)
check("Z3 up/charm charge difference", (2 - 0) % 3, 2, rel=1e-12)
check("bulk over charm sector-energy ratio", E_BULK / E_CHARM, 30.4, rel=3e-3)
check("EW sin^2 theta formula", 1 / (1 + PI), 0.2415, rel=2e-4)
check("EW identity Ee/sinW", E_EDGE / (1 / (1 + PI)), E_CHARM, rel=1e-12)
check("top mixing route value", E_BULK / E_TOP, 0.6455, rel=5e-5)
check("Route C exp(MU gap)=mu0", math.exp(MU * GAP), MU0, rel=1e-12)

trace2 = E_EDGE**2 + E_BOUNDARY**2 + E_BULK**2
check(
    "Route A Tr(Xsec^2)",
    trace2,
    162710.0,
    rel=1e-3,
    detail="Expected fail: pi^2+pi^4+16*pi^6 is about 15,489.506, not 162,710.",
)
check(
    "Route A 16*pi^6 component",
    16 * PI**6,
    162603.0,
    rel=1e-3,
    detail="Expected fail: 16*pi^6 is about 15,382.227, not 162,603.",
)
cube_root_trace = PI ** (1 / 3) + PI ** (2 / 3) + (4 * PI**3) ** (1 / 3)
check(
    "rational power trace Xsec^(1/3)",
    cube_root_trace,
    8.438,
    rel=1e-3,
    detail="Expected fail: direct evaluation gives about 8.5966.",
)

real_root = root_mu_equals_one()
check(
    "real k where mu_k=1",
    real_root,
    189.0,
    rel=1e-2,
    detail="Expected fail: solving the displayed mu_k formula gives k≈528.07.",
)
check(
    "mu_189 estimate",
    mu_k(189),
    0.64,
    rel=5e-2,
    detail="Expected fail: mu_189 is about 2.76; even the first term 16*pi^3/(189+4) is about 2.57.",
)

record(
    "moment-ratio proof covers all k,k' rather than mainly the k=0 reduction",
    False,
    computed="the written proof demonstrates the k=0 reduction mu_k'=1; it does not analyze arbitrary integer pairs k,k'",
    claimed="for all k,k' >= 0, mu_k/mu_k' = mu0 has no solution",
    detail="Expected proof-audit fail: the conclusion is plausible numerically but not proved by the displayed monotonic argument alone.",
)
record(
    "OP-B-Mass is explicitly left open",
    "Axiom~(MA) derivable from corpus" in TEX and "\\textbf{Open}" in TEX,
    computed="status table marks Axiom (MA) derivation open",
    claimed="bounded, not closed",
)

print(f"\n{'='*60}\nRESULT: {PASS} PASS / {FAIL} FAIL")
sys.exit(0 if FAIL == 0 else 1)