#!/usr/bin/env python3
"""
verify_P038.py -- Addendum 38: Quark sector and neutrino structure.

This verifier checks the explicit arithmetic in 38_Addendum_QuarkNeutrino.tex:
the density constants, self-lensing weights, quark spectral energies, and the
charm/top/Cabibbo/bottom numerical conjectures.  It also records one proof-audit
issue: Z3 cyclic symmetry alone does not uniquely determine the full TBM matrix.
"""

import math
import sys
from pathlib import Path

sys.path.insert(0, str(Path(__file__).resolve().parent))

PASS = FAIL = 0
_N = 0


def check(n, desc, cond):
    global PASS, FAIL
    ok = bool(cond)
    PASS += ok
    FAIL += not ok
    print(f"  [{'PASS' if ok else 'FAIL'}] {n:>2}. {desc}")
    return ok


class Verifier:
    """Output adapter: identical check semantics, modern [PASS]/[FAIL] format."""

    def __init__(self, name):
        print(name)

    def check(self, 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 self.record(label, ok, computed, claimed, detail, err_detail)

    def record(self, label, ok, computed="", claimed="", detail="", err_detail=""):
        global _N
        _N += 1
        desc = label + (f" -- {detail}" if detail else "")
        check(_N, desc, ok)
        if computed != "" or claimed != "":
            print(f"        computed: {computed}")
            print(f"        claimed : {claimed}")
        if err_detail:
            print(f"        {err_detail}")
        return ok

    def summary(self):
        print(f"\n{'='*60}\nRESULT: {PASS} PASS / {FAIL} FAIL")
        return 1 if FAIL else 0


v = Verifier("P038 -- Quark Sector and Neutrino Structure")

PI = math.pi
M_E = 0.51099895
MU0 = 4 * PI**3 + PI**2 + PI
MU1 = 16 * PI**3 / 5 + 3 * PI**2 / 4 + 2 * PI / 3
MU = MU1 / MU0
MU0_SUB = PI**2 + PI
MU1_SUB = 3 * PI**2 / 4 + 2 * PI / 3
E_LOW = 13.177
E_HIGH = 13.617

QUARK_MASSES = {
    "u": 2.16,
    "d": 4.67,
    "s": 93.4,
    "c": 1270.0,
    "b": 4180.0,
    "t": 172690.0,
}


def spectral_e(m_mev: float) -> float:
    return PI + math.log(m_mev / M_E) / MU


def mass_from_e(energy: float) -> float:
    return M_E * math.exp(MU * (energy - PI))


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


v.check("mu0 = 4*pi^3 + pi^2 + pi", MU0, 137.036, rel=3e-6)
v.check("mu1 = 16*pi^3/5 + 3*pi^2/4 + 2*pi/3", MU1, 108.717, rel=3e-6)
v.check("MU = mu1/mu0", MU, 0.7933, rel=1e-4)
v.check("mu0_sub = pi^2 + pi", MU0_SUB, 13.011, rel=2e-5)
v.check("mu1_sub = 3*pi^2/4 + 2*pi/3", MU1_SUB, 9.497, rel=5e-5)

lambda_tau = MU1 / (MU0 + MU1)
e_tau_from_endpoints = (1 - lambda_tau) * E_LOW + lambda_tau * E_HIGH
v.check("lambda_tau", lambda_tau, 0.4424, rel=5e-5)
v.check(
    "E_self(lambda_tau) from stated endpoints",
    e_tau_from_endpoints,
    13.419,
    rel=5e-4,
    detail="Expected fail: endpoints 13.177 and 13.617 with lambda_tau=0.4424 give about 13.372; 13.419 would require lambda=0.55.",
)

lambda_sub = MU1_SUB / (MU0_SUB + MU1_SUB)
e_sub = (1 - lambda_sub) * E_LOW + lambda_sub * E_HIGH
sub_ratio = math.exp(MU * (e_sub - PI))
v.check("lambda_sub", lambda_sub, 0.422, rel=2e-4)
v.check("E_self(lambda_sub)", e_sub, 13.36, rel=3e-4)
v.check("sub-fold predicted mass ratio", sub_ratio, 3300, rel=1e-2)
v.check("sub-fold predicted mass in MeV", M_E * sub_ratio, 1690, rel=1e-2)
v.record("muon energy below stated self-lensing range", PI**2 < E_LOW, f"{PI**2:.8f}", "< 13.177")

lambda_mu = PI / (1 + PI)
emu_from_ratio = PI + math.log(206.768) / MU
v.check("lambda_mu = pi/(1+pi)", lambda_mu, 0.759, rel=7e-4)
v.check("muon mass formula ratio from pi^2", math.exp(MU * (PI**2 - PI)), 208.0, rel=5e-4)
v.check("muon mass residual percent", residual(math.exp(MU * (PI**2 - PI)), 206.768), 0.60, rel=1e-2)
v.check("Hopf-holonomy delta E", emu_from_ratio - PI**2, -0.007584, rel=5e-4)
e_ct = (PI + PI**2) / 2
v.check("Clifford torus energy", e_ct, 6.506, rel=7e-5)
v.check("Clifford torus mass ratio", math.exp(MU * (e_ct - PI)), 14.4, rel=2e-3)
v.check("Clifford torus mass in MeV", mass_from_e(e_ct), 7.4, rel=5e-3)

v.check("TBM sin^2 theta12", 1 / 3, 0.333, rel=2e-3)
v.check("TBM sin^2 theta23", 1 / 2, 0.500, rel=1e-12)
v.check("TBM theta13", 0.0, 0.0, abs_tol=0.0)
v.check("theta12 TBM absolute residual percent", abs(residual(1 / 3, 0.307)), 8.0, rel=8e-2)
v.check("theta23 TBM absolute residual percent", abs(residual(1 / 2, 0.545)), 9.0, rel=1e-1)
v.record(
    "Z3 cyclic symmetry alone uniquely fixes TBM",
    False,
    "Z3 fixes the trimaximal eigenvector but leaves a 2D subspace",
    "unique TBM",
    "Expected proof-audit fail: an extra Z2/S4-style constraint is needed for the full TBM basis.",
)

e_nu = PI + math.log(0.05 / (M_E * 1_000_000)) / MU
v.check("representative 0.05 eV neutrino energy", e_nu, -17.2, rel=5e-3)
v.record("neutrino representative energy lies below electron floor", e_nu < PI, f"{e_nu:.6f}", "< pi")

claimed_quark_e = {
    "u": 4.959,
    "d": 5.931,
    "s": 9.701,
    "c": 12.996,
    "b": 14.498,
    "t": 19.187,
}
for quark, claim in claimed_quark_e.items():
    tol = 3e-4 if quark == "s" else 1e-4
    v.check(
        f"E_{quark} from PDG mass",
        spectral_e(QUARK_MASSES[quark]),
        claim,
        rel=tol,
        detail="Expected fail: with the paper's constants, E_s is about 9.7066, not 9.701." if quark == "s" else "",
    )

ed_log = math.log(MU1) / MU
ed_pdg = spectral_e(QUARK_MASSES["d"])
v.check("down conjecture E_d = ln(mu1)/MU", ed_log, 5.910, rel=3e-5)
v.check("down conjecture absolute residual percent", abs(residual(ed_log, ed_pdg)), 0.34, rel=2e-2)
v.check("ln(mu0)/ln(mu1)", math.log(MU0) / math.log(MU1), 1.049, rel=5e-4)

m_charm_pred = mass_from_e(MU0_SUB)
v.check("charm energy target pi^2 + pi", MU0_SUB, 13.011, rel=2e-5)
v.check("charm predicted mass", m_charm_pred, 1285, rel=5e-4)
v.check("charm mass residual percent", residual(m_charm_pred, QUARK_MASSES["c"]), 1.2, rel=2e-2)

top_charm_ratio = QUARK_MASSES["t"] / QUARK_MASSES["c"]
top_charm_gap_pred = math.log(MU0) / MU
top_charm_gap_pdg = spectral_e(QUARK_MASSES["t"]) - spectral_e(QUARK_MASSES["c"])
v.check("top/charm mass ratio", top_charm_ratio, 135.98, rel=5e-5)
v.check("top/charm ratio residual percent", residual(top_charm_ratio, MU0), -0.77, rel=1e-2)
v.check("top/charm predicted energy gap", top_charm_gap_pred, 6.202, rel=3e-4)
v.check("top/charm PDG energy gap", top_charm_gap_pdg, 6.191, rel=3e-4)
v.check("top/charm gap difference", top_charm_gap_pred - top_charm_gap_pdg, 0.011, rel=2e-1)

cabibbo = math.sin(PI / 14)
v.check("Cabibbo sin(pi/14)", cabibbo, 0.22252, rel=5e-6)
v.check("Cabibbo residual percent", residual(cabibbo, 0.22431), -0.80, rel=1e-2)
v.check("naive Vcb = lambda^2", cabibbo**2, 0.0495, rel=5e-4)
v.check("naive Vub = lambda^3", cabibbo**3, 0.0110, rel=2e-3)

e_bottom = spectral_e(QUARK_MASSES["b"])
e_bottom_fano = PI ** (7 / 3)
v.check("bottom Fano energy pi^(7/3)", e_bottom_fano, 14.4549, rel=5e-6)
v.check("bottom empirical energy", e_bottom, 14.4980, rel=1e-5)
v.check("bottom energy residual percent", residual(e_bottom_fano, e_bottom), -0.30, rel=2e-2)

fermion_energies = {
    "e": PI,
    "mu": emu_from_ratio,
    "tau": spectral_e(1776.86),
    **{q: spectral_e(m) for q, m in QUARK_MASSES.items()},
}
lattice_hits = []
lattice_misses = []
for n in range(0, 12):
    value = PI ** (n / 3)
    best_name, best_energy = min(
        fermion_energies.items(), key=lambda item: abs(value / item[1] - 1)
    )
    miss = abs(residual(value, best_energy))
    if n in {3, 6, 7}:
        lattice_hits.append((n, miss, best_name))
    else:
        lattice_misses.append((n, miss, best_name))
v.record("pi^(n/3) claimed hits are within 1 percent", all(miss < 1 for _, miss, _ in lattice_hits), lattice_hits, "n=3,6,7")
v.record("other pi^(n/3) entries miss by more than 7 percent through n=11", all(miss > 7 for _, miss, _ in lattice_misses), lattice_misses, ">7%")

sys.exit(v.summary())