#!/usr/bin/env python3
"""
verify_P058.py -- Addendum 58: seesaw scale closure status.

This verifier checks the numerical backbone of P58:
  * TOE constants and spectral coordinates for v, m_nu3, and M_R
  * the Phase 5a mass-squared ratio r
  * the anchor-equivalence formulas for masses and Sigma m_nu
  * the natural-expression search tables
  * the Option B exclusion arithmetic

Several failures are intentional audit findings: P58 often makes the correct
status claim that M_R remains open, but it mixes mass ratios with
mass-squared ratios in a few formulas and contains some stale table arithmetic.
"""

import math
import sys
from pathlib import Path

sys.path.insert(0, str(Path(__file__).resolve().parent))
from verify_common import CheckResult, Verifier


class ModernVerifier(Verifier):
    """Local adapter: tolerance logic inherited byte-identical from
    verify_common.Verifier; only the output layer is modernised."""

    def __init__(self, name: str) -> None:
        self.name = name
        self.results = []
        print(name)

    def record(self, label, ok, computed="", claimed="", detail=""):
        self.results.append(CheckResult(label, ok, computed, claimed, detail))
        n = len(self.results)
        desc = f"{label} -- {detail}" if detail else label
        print(f"  [{'PASS' if ok else 'FAIL'}] {n:>2}. {desc}")
        if computed != "" or claimed != "":
            print(f"        computed: {computed}")
            print(f"        claimed : {claimed}")
        return ok

    def summary(self):
        passed = sum(r.ok for r in self.results)
        failed = len(self.results) - passed
        print(f"\n{'='*60}\nRESULT: {passed} PASS / {failed} FAIL")
        return 1 if failed else 0


v = ModernVerifier("P058 -- Seesaw Scale Closure After Phase 5a")

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

PI = math.pi
M_E_MEV = 0.51100
V_EW_GEV = 246.22
MNU3_MEV = 49.5e-9
MNU3_GEV = MNU3_MEV / 1000
MR_TABLE_GEV = 1.24e15
DM31_EV2 = 2.453e-3
DM31_UNCERT_EV2 = 0.034e-3

MU0 = 4 * PI**3 + PI**2 + PI
MU1 = 16 * PI**3 / 5 + 3 * PI**2 / 4 + 2 * PI / 3
MU = MU1 / MU0

LAMBDA = math.sin(PI / 14)
Y1 = LAMBDA
Y2 = math.sin(2 * PI / 14)


def spectral_from_gev(mass_gev: float) -> float:
    return PI + math.log((mass_gev * 1000) / M_E_MEV) / MU


def spectral_from_mev(mass_mev: float) -> float:
    return PI + math.log(mass_mev / M_E_MEV) / MU


def mass_gev_from_e(e: float) -> float:
    return (M_E_MEV * math.exp(MU * (e - PI))) / 1000


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


v.check("mu0 = 4*pi^3 + pi^2 + pi", MU0, 137.036, rel=5e-6)
v.check("mu1 continuous moment", MU1, 108.717, rel=5e-6)
v.check("MU = mu1/mu0", MU, 0.79334, rel=5e-5)

e_v = spectral_from_gev(V_EW_GEV)
e_nu3_from_49_5_mev = spectral_from_mev(MNU3_MEV)
e_r_from_table_mr = spectral_from_gev(MR_TABLE_GEV)
mr_from_49_5 = V_EW_GEV**2 / MNU3_GEV
e_r_from_seesaw_49_5 = spectral_from_gev(mr_from_49_5)
reflection_from_49_5 = 2 * e_v - e_nu3_from_49_5_mev

v.check("E_v from v = 246.22 GeV", e_v, 19.636, abs_tol=0.002)
v.check(
    "E_nu3 from listed m_nu3 = 49.5 meV",
    e_nu3_from_49_5_mev,
    -17.231,
    abs_tol=0.005,
    detail="Expected fail: 49.5 meV gives about -17.215; -17.231 corresponds to a slightly lower mass or to the rounded reflection chain.",
)
v.check("E_R from listed M_R = 1.24e15 GeV", e_r_from_table_mr, 56.502, abs_tol=0.005)
v.check(
    "M_R from listed v and m_nu3",
    mr_from_49_5,
    1.24e15,
    rel=5e-3,
    detail="Expected fail: v^2/(49.5 meV) gives about 1.224e15 GeV, not 1.24e15.",
)
v.check(
    "spectral reflection using listed v and m_nu3",
    reflection_from_49_5,
    56.503,
    abs_tol=0.005,
    detail="Expected fail: using the listed 49.5 meV input gives about 56.486, not 56.503.",
)
v.check("E_R from seesaw M_R using listed m_nu3", e_r_from_seesaw_49_5, reflection_from_49_5, abs_tol=1e-12)

r_toe = (Y2**4 - Y1**4) / (1 - Y1**4)
v.check("y1 = sin(pi/14)", Y1, 0.22252, rel=5e-6)
v.check("y2 = sin(2pi/14)", Y2, 0.43388, rel=1e-5)
v.check("y1^4", Y1**4, 0.002455, rel=2e-3)
v.check("y2^4", Y2**4, 0.035477, rel=2e-3)
v.check("r_TOE mass-squared ratio", r_toe, 0.03307, rel=5e-4)
v.check("r_TOE vs r_PDG relative excess", pct(r_toe, 0.03070), 7.7, rel=5e-3)

v.record(
    "leading mass spectrum should use sin^2, not sin^4",
    "m_{\\nu_2}^{(0)} &= \\sin^4(2\\pi/14)" not in TEX,
    "TeX prints sin^4 for m_nu2 and m_nu1",
    "m_nu_k = y_k^2 v^2/M_R",
    "Expected fail: the paper's own earlier seesaw lemma gives masses proportional to y_k^2; sin^4 enters only after squaring masses for Delta m^2.",
)
v.check(
    "m_nu1/m_nu3 mass ratio",
    Y1**2,
    Y1**4,
    rel=1e-3,
    detail="Expected fail: printed sin^4 mass coefficient is the squared-mass coefficient.",
)
v.check(
    "m_nu2/m_nu3 mass ratio",
    Y2**2,
    Y2**4,
    rel=1e-3,
    detail="Expected fail: printed sin^4 mass coefficient is the squared-mass coefficient.",
)
v.record(
    "Sigma m_nu proof uses the same powers as the proposition",
    "m_{\\nu_3}(1 + \\sin^4(\\pi/14) + \\sin^4(2\\pi/14))" not in TEX,
    "proof uses sin^4 factor",
    "proposition states 1 + y1^2 + y2^2",
    "Expected fail: Sigma m_nu should use mass ratios y_k^2, not squared-mass ratios y_k^4.",
)
v.check("Sigma mass factor using sin^2", 1 + Y1**2 + Y2**2, 1.23777, rel=5e-5)
v.check(
    "Sigma mass factor if proof's sin^4 is used",
    1 + Y1**4 + Y2**4,
    1.23777,
    rel=1e-3,
    detail="Expected fail: the proof's sin^4 factor gives about 1.03789.",
)

pdg_r = 0.03070
delta_r = pdg_r - r_toe
dr_dy2 = 4 * Y2**3 / (1 - Y1**4)
delta_y2 = delta_r / dr_dy2
delta_y2_frac = delta_y2 / Y2
c_peirce_required = -delta_y2_frac / Y1**2
v.check("required delta r / r_TOE", 100 * delta_r / r_toe, -7.16, rel=5e-3)
v.check(
    "required delta y2/y2",
    100 * delta_y2_frac,
    3.5,
    rel=5e-2,
    detail="Expected fail: the linearized correction is about -1.67%, not +3.5%.",
)
v.record("required c_Peirce lies in stated [0.3,0.4]", 0.3 <= c_peirce_required <= 0.4, f"{c_peirce_required:.6f}", "[0.3,0.4]")

target_abs_enu = 17.231
search_candidates = {
    "2*pi^2/MU": 2 * PI**2 / MU,
    "pi^2 + pi": PI**2 + PI,
    "pi^2 + 2*pi": PI**2 + 2 * PI,
    "pi^2 + 2*pi + 1/(2*pi)": PI**2 + 2 * PI + 1 / (2 * PI),
    "mu1*pi/(pi^2 + 3*pi)": MU1 * PI / (PI**2 + 3 * PI),
    "mu1/(2*pi^2)": MU1 / (2 * PI**2),
    "pi^2 + 2*pi + 1/pi": PI**2 + 2 * PI + 1 / PI,
    "pi^2 + 2*pi + pi*MU/(pi+1)": PI**2 + 2 * PI + PI * MU / (PI + 1),
}
v.check("search table: 2*pi^2/MU", search_candidates["2*pi^2/MU"], 24.71, rel=1e-2)
v.check(
    "search table: mu1*pi/(pi^2 + 3*pi)",
    search_candidates["mu1*pi/(pi^2 + 3*pi)"],
    19.28,
    rel=5e-3,
    detail="Expected fail: the displayed formula gives about 17.70.",
)
best_search_error = min(abs(value / target_abs_enu - 1) for value in search_candidates.values())
v.record(
    "search table supports 'no expression within 6%'",
    best_search_error > 0.06,
    f"best listed relative error {100 * best_search_error:.2f}%",
    ">6%",
    "Expected fail: the paper's own table includes candidates at about 5.3%, 4.4%, and 2.8%.",
)
v.record(
    "pi^2 + 2*pi is the closest listed simple expression",
    abs(search_candidates["pi^2 + 2*pi"] / target_abs_enu - 1) <= best_search_error + 1e-12,
    f"pi^2+2pi error {100 * abs(search_candidates['pi^2 + 2*pi'] / target_abs_enu - 1):.2f}%; best listed {100 * best_search_error:.2f}%",
    "closest",
    "Expected fail: later rows in the same table are closer.",
)

e_tilde_from_claimed = 2 * PI + target_abs_enu
e_tilde_from_listed_mnu = 2 * PI - e_nu3_from_49_5_mev
v.check("mirror E_tilde from claimed |E_nu3|", e_tilde_from_claimed, 23.513, abs_tol=0.002)
v.check("mirror E_tilde from listed 49.5 meV", e_tilde_from_listed_mnu, 23.513, abs_tol=0.005)
v.check(
    "mirror table: pi^2 + pi^2/4",
    PI**2 + PI**2 / 4,
    14.81,
    rel=5e-3,
    detail="Expected fail: pi^2 + pi^2/4 is about 12.34; the row's arithmetic double-counts.",
)
mirror_candidates = [
    PI**2 + PI**2 / 4,
    3 * PI**2 / 4 + 2 * PI,
    MU1 / (3 * PI**2),
    PI**2 + PI + PI / MU,
    math.sqrt(MU0),
    3 * PI,
    3 * PI + PI**2,
    4 * PI - 1,
    3 * PI + PI**2 / 2,
    4 * PI + PI**2 / 4,
]
best_mirror_error = min(abs(value / e_tilde_from_claimed - 1) for value in mirror_candidates)
v.record("no mirror candidate within 10%", best_mirror_error > 0.10, f"best relative error {100 * best_mirror_error:.2f}%", ">10%")
v.record(
    "mirror relation gives positive |E|",
    2 * PI - e_tilde_from_claimed > 0,
    f"2*pi - E_tilde = {2 * PI - e_tilde_from_claimed:.6f}",
    "|E_nu3|",
    "Expected fail: if E_tilde = 2*pi - E_nu3, then |E_nu3| = E_tilde - 2*pi, not 2*pi - E_tilde.",
)

e_b = MU1 / 2 + PI * MU
mr_b = mass_gev_from_e(e_b)
m3_b_mev = V_EW_GEV**2 / mr_b * 1e12
m1_b_sin2 = Y1**2 * m3_b_mev
m2_b_sin2 = Y2**2 * m3_b_mev
m1_b_sin4 = Y1**4 * m3_b_mev
m2_b_sin4 = Y2**4 * m3_b_mev
dm31_b_ev2 = (m3_b_mev / 1000) ** 2 - (m1_b_sin2 / 1000) ** 2
v.check("Option B energy mu1/2 + pi*MU", e_b, 56.85, rel=5e-4)
v.check("Option B M_R", mr_b, 1.64e15, rel=5e-3)
v.check("Option B m_nu3", m3_b_mev, 37.0, rel=5e-3)
v.check(
    "Option B m_nu1 using printed sin^4 formula",
    m1_b_sin4,
    0.91,
    rel=2e-2,
    detail="Expected fail: sin^4 gives about 0.091 meV; correct sin^2 gives about 1.84 meV.",
)
v.check(
    "Option B m_nu2 using printed sin^4 formula",
    m2_b_sin4,
    6.89,
    rel=2e-2,
    detail="Expected fail: sin^4 gives about 1.31 meV; the quoted 6.89 meV uses sin^2.",
)
v.check("Option B m_nu2 if corrected sin^2 is used", m2_b_sin2, 6.89, rel=2e-2)
v.check("Option B Delta m31", dm31_b_ev2, 1.37e-3, rel=5e-3)
v.check("Option B Delta m31 discrepancy", pct(dm31_b_ev2, DM31_EV2), -44.0, rel=1e-2)

sqrt_dm31_gev = math.sqrt(DM31_EV2) * 1e-9
mr_from_dm31 = V_EW_GEV**2 / sqrt_dm31_gev
mr_uncert = mr_from_dm31 * 0.5 * DM31_UNCERT_EV2 / DM31_EV2
e_uncert = (0.5 * DM31_UNCERT_EV2 / DM31_EV2) / MU
v.check("M_R from displayed Delta m31 shortcut", mr_from_dm31, 1.24e15, rel=2e-2)
v.check(
    "propagated M_R uncertainty from Delta m31",
    mr_uncert / 1e15,
    0.02,
    rel=1e-1,
    detail="Expected fail: the displayed Delta m31 uncertainty gives about +/-0.0085e15 GeV, not +/-0.02e15.",
)
v.check("propagated E_R absolute uncertainty", e_uncert, 0.0088, rel=1e-2)

v.record(
    "quartic section does not reuse trilinear T as a 4-linear form",
    "N_4(X) =\nT(X,X,X,X)" not in TEX and "N_4(X) = T(X,X,X,X)" not in TEX,
    "TeX defines N_4(X) = T(X,X,X,X)",
    "T is trilinear earlier in the corpus",
    "Expected fail: the fourth-order construction needs its own 4-linear invariant; writing the cubic trilinear T with four arguments is not a valid definition.",
)

sys.exit(v.summary())