#!/usr/bin/env python3
"""
verify_P040.py -- Addendum 40: J3(O) spectral table and PMNS corrections.

This verifier checks the spectral-table masses, PMNS Phase 3.5 numerics, and
Higgs-VEV candidate arithmetic in 40_Addendum_SpectralTable.tex.  The mass
table and PMNS formulas mostly reproduce.  The flagged issues are bookkeeping
and internal arithmetic problems: the derivation-status slot counts do not
match the listed slots, the VEV residual sign is reversed under the table's
own convention, and the tree-level Weinberg VEV route uses the wrong square
root factor.
"""

from __future__ import annotations

import math
import re
import sys
from pathlib import Path

PASS = FAIL = 0
_N = 0


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


class Verifier:
    """Same check semantics as verify_common.Verifier; modern output style."""

    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 = abs(computed - claimed)
            err_detail = f"abs err={err:.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, err_detail + (f"; {detail}" if detail else ""))

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

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


v = Verifier("P040 -- Spectral Table and PMNS Corrections")

ROOT = Path(__file__).resolve().parents[1]
TEX = (ROOT / "40_Addendum_SpectralTable.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
ME_MEV = 0.511
MZ = 91.188
V_PDG = 246.22


def mass_mev(energy: float) -> float:
    return ME_MEV * math.exp(MU * (energy - PI))


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


theta_c = PI / 14
sin2_c = math.sin(theta_c) ** 2

v.check("mu0", MU0, 137.036, rel=3e-6)
v.check("mu1", MU1, 108.717, rel=5e-6)
v.check("MU", MU, 0.7933, rel=6e-5)
v.check("sin(theta_C)", math.sin(theta_c), 0.22252, rel=5e-5)

v.check("electron mass at E=pi", mass_mev(PI), 0.511, rel=1e-12)
v.check("muon mass at E=pi^2", mass_mev(PI**2), 106.3, rel=5e-4)
v.check("muon residual", pct(mass_mev(PI**2), 105.7), 0.60, rel=8e-2)
down_energy = math.log(MU1) / MU
v.check("down/Cabibbo-direction mass", mass_mev(down_energy), 4.60, rel=2e-3)
v.check("down residual", pct(mass_mev(down_energy), 4.67), -1.6, rel=2e-2)
charm_energy = PI**2 + PI
v.check("charm mass", mass_mev(charm_energy), 1285, rel=5e-4)
v.check("charm residual", pct(mass_mev(charm_energy), 1270), 1.2, rel=2e-2)
top_energy = PI**2 + PI + math.log(MU0) / MU
v.check("top mass", mass_mev(top_energy) / 1000, 176, rel=2e-3)
v.check("top residual vs 173 GeV", pct(mass_mev(top_energy) / 1000, 173), 1.9, rel=7e-2)
bottom_energy = PI ** (7 / 3)
v.check("bottom mass", mass_mev(bottom_energy), 4040, rel=5e-4)
v.check("bottom residual", pct(mass_mev(bottom_energy), 4180), -3.3, rel=2e-2)

slot_numbers = [int(n) for n in re.findall(r"^\s*(\d+)\s*&", TEX, flags=re.MULTILINE)]
v.check("spectral table slot count", len(slot_numbers), 27, rel=0)
known_formula_slots = [1, 2, 3, 5, 9, 10, 11, 13, 14, 15, 16, 22, 23, 24, 25, 26, 27]
candidate_formula_slots = [2, 20, 25, 26, 27]
open_slots_listed = [4, 6, 7, 8, 12, 17, 18, 19, 21]
v.record(
    "derivation-status known-formula slot count is 16",
    len(known_formula_slots) == 16,
    computed=f"listed formula slots={len(known_formula_slots)}",
    claimed="16 carry known TOE energy formulas",
    detail="Expected bookkeeping fail: the slot ranges named in the sentence contain 17 slots.",
)
v.record(
    "derivation-status candidate-formula slot count is 4",
    len(candidate_formula_slots) == 4,
    computed=f"slot 2 + slot 20 + slots 25-27 = {len(candidate_formula_slots)} slots",
    claimed="4 carry candidate formulas",
    detail="Expected bookkeeping fail if counted by slots rather than formula classes.",
)
v.record(
    "derivation-status open slot count is 7",
    len(open_slots_listed) == 7,
    computed=f"listed open slots={len(open_slots_listed)}",
    claimed="7 remain open",
    detail="Expected bookkeeping fail: slots 4, 6-8, 12, 17-19, and 21 total 9 slots.",
)

theta13 = (4 / 9) * sin2_c
theta23 = 15 / 28
theta12 = math.cos(theta_c) ** 2 / 3
v.check("PMNS theta13 conjecture", theta13, 0.022007, rel=5e-5)
v.check("PMNS theta13 residual", pct(theta13, 0.02195), 0.26, rel=5e-2)
v.check("PMNS theta23 conjecture", theta23, 0.53571, rel=1e-5)
v.check("PMNS theta23 residual", pct(theta23, 0.546), -1.88, rel=8e-3)
v.check("PMNS theta12 conjecture", theta12, 0.31683, rel=1e-5)
v.check("PMNS theta12 residual", pct(theta12, 0.307), 3.20, rel=5e-3)
v.check("TBM theta12 residual", pct(1 / 3, 0.307), 8.6, rel=5e-3)
v.check("TBM theta23 residual", pct(1 / 2, 0.546), -8.4, rel=5e-3)
v.check("TBM theta13 residual", pct(0.0, 0.02195), -100.0, rel=1e-12)

v_vev_candidate = math.e * MZ
v.check("Higgs VEV candidate e*MZ", v_vev_candidate, 247.87, rel=2e-5)
v.check(
    "Higgs VEV candidate residual sign",
    pct(v_vev_candidate, V_PDG),
    -0.67,
    rel=2e-2,
    detail="Expected fail: with calc-minus-PDG convention the residual is +0.672%, not -0.67%.",
)
ev_energy = PI + math.log((V_PDG * 1000) / ME_MEV) / MU
v.check("Higgs VEV energy", ev_energy, 19.636, rel=3e-5)
v.record(
    "slot-12 to slot-20 one-block shift is code-verifiable in P40",
    False,
    computed="slot 12 has no energy formula in the table",
    claimed="E_v = E(slot 12) + 1/MU ≈ 19.644",
    detail="Expected status fail: this numerical interpretation uses an unspecified slot-12 energy.",
)

tree_factor = math.sqrt((1 + PI) / PI)
tree_vev = MZ * tree_factor
v.check(
    "tree Weinberg VEV square-root factor",
    tree_factor,
    2.032,
    rel=2e-3,
    detail="Expected fail: sqrt((1+pi)/pi)≈1.148, not 2.032.",
)
v.check(
    "tree Weinberg VEV value",
    tree_vev,
    185.3,
    rel=2e-3,
    detail="Expected fail: MZ*sqrt((1+pi)/pi)≈104.7 GeV, not 185.3 GeV.",
)

e_u_from_mass = PI + math.log(2.16 / ME_MEV) / MU
v.check("u-quark energy from 2.16 MeV", e_u_from_mass, 4.959, rel=1e-4)
v.check("u candidate pi^2/2 energy residual", pct(PI**2 / 2, e_u_from_mass), -0.49, rel=3e-2)
e_s_from_mass = PI + math.log(93.4 / ME_MEV) / MU
v.check("s-quark energy from 93.4 MeV", e_s_from_mass, 9.707, rel=5e-5)

sys.exit(v.summary())