#!/usr/bin/env python3
"""
verify_P083.py -- Addendum 83: Light Quark Masses.

This verifier checks the light-quark spectral-energy and mass claims in
83_Addendum_LightQuarkMasses.tex.

It intentionally flags the internal inconsistency in Theorem D: the abstract
states that the direct TOE formulas give E_d - E_u ~= 0.9774 with an error
below 0.002 relative to 2 lambda^2 pi^2, but the paper's own formulas give
E_d - E_u ~= 0.9441. The body later acknowledges this 3-4% discrepancy, so the
problem is the abstract/theorem statement, not the u, d, s mass arithmetic.
"""

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("P083 -- Light Quark Masses")

PI = math.pi
M_E = 0.511
M_MUON = 105.66
MU0 = 4 * PI**3 + PI**2 + PI
MU1 = 16 * PI**3 / 5 + 3 * PI**2 / 4 + 2 * PI / 3
MU = MU1 / MU0
LAM = math.sin(PI / 14)
N_IM = 7
N_C = 3

PDG = {
    "u": 2.16,
    "d": 4.67,
    "s": 93.4,
}

E = {
    "u": PI ** (N_IM / (N_IM - N_C + 1)),
    "d": math.log(MU1) / MU,
    "s": PI**2 - 1 / N_IM,
}


def mass_from_e(e):
    return M_E * math.exp(MU * (e - PI))


def e_from_mass(m):
    return PI + math.log(m / M_E) / MU


def residual_pct(pred, actual):
    return 100 * (pred - actual) / actual


check("mu0 = 4*pi^3 + pi^2 + pi", MU0, 137.036, rel=3e-6)
check("mu1 = 16*pi^3/5 + 3*pi^2/4 + 2*pi/3", MU1, 108.717, rel=3e-6)
check("MU = mu1/mu0", MU, 0.7933, rel=1e-4)
check("lambda = sin(pi/14)", LAM, 0.22252, rel=5e-6)

claimed_empirical_e = {
    "u": 4.9582,
    "d": 5.9309,
    "s": 9.7073,
}

for q, claimed in claimed_empirical_e.items():
    check(f"E_{q} empirical from PDG mass", e_from_mass(PDG[q]), claimed, rel=2e-4)

claimed_spectral_e = {
    "u": 4.9672,
    "d": 5.9114,
    "s": 9.7267,
}

for q, claimed in claimed_spectral_e.items():
    check(f"E_{q} TOE spectral formula", E[q], claimed, rel=4e-4)

pred_mass = {q: mass_from_e(E[q]) for q in E}
claimed_masses = {
    "u": 2.19,
    "d": 4.59,
    "s": 95.1,
}

for q, claimed in claimed_masses.items():
    check(f"m_{q} from spectral formula", pred_mass[q], claimed, rel=1.5e-2)

actual_residuals = {q: residual_pct(pred_mass[q], PDG[q]) for q in E}
for q, residual in actual_residuals.items():
    record(
        f"{q} residual is within the paper's stated 2% band",
        abs(residual) < 2.0,
        f"{residual:+.4f}%",
        "<2%",
    )

direct_gap = E["d"] - E["u"]
weyl_gap = 2 * LAM**2 * PI**2
pdg_gap = e_from_mass(PDG["d"]) - e_from_mass(PDG["u"])
epsilon = direct_gap - weyl_gap
gap_discrepancy_pct = 100 * (direct_gap / weyl_gap - 1)

check("direct E_d - E_u gap from TOE formulas", direct_gap, 0.9442, rel=5e-4)
check("Weyl-step gap 2*lambda^2*pi^2", weyl_gap, 0.9783, rel=1.5e-3)
check("PDG-derived E_d - E_u gap", pdg_gap, 0.9727, rel=1.5e-3)
check("Weyl-step gap vs PDG residual percent", 100 * (weyl_gap / pdg_gap - 1), 0.58, rel=2e-1)
check(
    "body discrepancy between direct gap and Weyl-step gap",
    abs(gap_discrepancy_pct),
    3.6,
    rel=7e-2,
)

check(
    "abstract/Theorem D direct-gap value",
    direct_gap,
    0.9774,
    rel=2e-3,
    detail="Expected fail: direct formulas give about 0.9441; 0.9774 is the Weyl-step value.",
)
record(
    "abstract/Theorem D epsilon bound |epsilon| < 0.002",
    abs(epsilon) < 0.002,
    f"{abs(epsilon):.8f}",
    "<0.002",
    "Expected fail: epsilon is direct_gap - 2*lambda^2*pi^2.",
)

md_mu_weyl = math.exp(MU * weyl_gap)
check("m_d/m_u Weyl-step formula", md_mu_weyl, 2.173, rel=2e-3)
check("PDG m_d/m_u ratio", PDG["d"] / PDG["u"], 2.162, rel=2e-3)
check("m_s/m_mu = exp(-MU/7)", math.exp(-MU / 7), 0.8927, rel=3e-4)
check("m_s from muon ratio", M_MUON * math.exp(-MU / 7), 94.3, rel=1e-3)
check("PDG m_s/m_mu ratio", PDG["s"] / M_MUON, 0.8840, rel=5e-4)
check("m_s/m_d formula ratio", math.exp(MU * (E["s"] - E["d"])), 20.61, rel=3e-3)
check("PDG m_s/m_d ratio", PDG["s"] / PDG["d"], 20.00, rel=5e-4)

record("light-quark spectral ordering", E["u"] < E["d"] < E["s"], f"{E}", "u < d < s")

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