#!/usr/bin/env python3
"""
verify_P061.py -- Addendum 61: Planck-Seesaw Spectral Gap.

This verifier checks the numerical backbone of
61_Addendum_PlanckSeesawGap.tex: Planck spectral coordinate, seesaw scale,
the Planck-seesaw gap, the exact E6 Coxeter conjecture inversion, and the
simple epsilon-candidate table.

It intentionally flags several internal consistency issues:
  * The theorem statement for E_R contains a MeV/GeV conversion error
    and gives E_R ~= 65.189. A later remark corrects this, but the theorem
    itself remains false.
  * The paper presents a "full precision" E_Pl/gap while adopting rounded
    legacy values; recomputing from its stated inputs shifts the numbers by
    more than the quoted tight uncertainty.
  * The abstract says +/-0.004 for the gap uncertainty, while the paper's
    own propagation gives about +/-0.0088.
  * The reduced Planck mass table uses an incorrect spectral coordinate.
"""

import math
import sys
from pathlib import Path

sys.path.insert(0, str(Path(__file__).resolve().parent))
class Verifier:
    """Output shim: identical tolerance semantics to verify_common.Verifier,
    modern [PASS]/[FAIL] check-line output format."""

    def __init__(self, name):
        self.PASS = 0
        self.FAIL = 0
        self.n = 0
        print(name)

    def _mark(self, ok, desc):
        self.n += 1
        if ok:
            self.PASS += 1
        else:
            self.FAIL += 1
        print(f"  [{'PASS' if ok else 'FAIL'}] {self.n:>2}. {desc}")

    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}%"
        self._mark(ok, label + (f" -- {detail}" if detail else ""))
        print(f"        computed: {computed}   claimed: {claimed}   ({err_detail})")
        return ok

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

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


v = Verifier("P061 -- Planck-Seesaw Spectral Gap")

PI = math.pi
M_E_MEV = 0.51100
M_PL_GEV = 1.22090e19
M_PL_REDUCED_GEV = 2.43532e18
V_EW_GEV = 246.22
DM31_EV2 = 2.453e-3
DM31_FRAC_UNCERT = 0.014
H_E6 = 12.0

MU0 = 4 * PI**3 + PI**2 + PI
MU1 = 16 * PI**3 / 5 + 3 * PI**2 / 4 + 2 * PI / 3
MU = MU1 / MU0
THETA_C = PI / 14
LAMBDA = math.sin(PI / 14)


def spectral_e_from_gev(mass_gev):
    return PI + math.log((mass_gev * 1000) / M_E_MEV) / MU


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


sqrt_dm31_ev = math.sqrt(DM31_EV2)
sqrt_dm31_gev = sqrt_dm31_ev * 1e-9
mr_from_listed_inputs_gev = V_EW_GEV**2 / sqrt_dm31_gev
mr_table_gev = 1.2425e15

e_pl = spectral_e_from_gev(M_PL_GEV)
e_pl_reduced = spectral_e_from_gev(M_PL_REDUCED_GEV)
e_r_from_listed_inputs = spectral_e_from_gev(mr_from_listed_inputs_gev)
e_r_from_table_mr = spectral_e_from_gev(mr_table_gev)
e_r_canonical = 56.502


v.check("mu0 = 4*pi^3 + pi^2 + pi", MU0, 137.036304, rel=5e-9)
v.check("mu1 = 16*pi^3/5 + 3*pi^2/4 + 2*pi/3", MU1, 108.716684, rel=5e-9)
v.check("MU = mu1/mu0", MU, 0.793338, rel=6e-6)

ratio_pl = (M_PL_GEV * 1000) / M_E_MEV
ln_pl = math.log(ratio_pl)
v.check("M_Pl/m_e ratio", ratio_pl, 2.38923e22, rel=5e-6)
v.check("ln(M_Pl/m_e)", ln_pl, 51.52794, abs_tol=1.5e-4)
v.check(
    "E_Pl from CODATA Planck mass and exact MU",
    e_pl,
    68.0946,
    abs_tol=1.5e-3,
    detail="Expected fail: exact arithmetic gives about 68.0919, not 68.0946/68.096.",
)

v.check("sqrt(Delta m31^2)", sqrt_dm31_ev, 0.04953, rel=1e-4)
v.check(
    "M_R from listed v_EW and Delta m31 inputs",
    mr_from_listed_inputs_gev,
    1.2425e15,
    rel=5e-3,
    detail="Expected fail: listed inputs give about 1.2240e15 GeV; 1.2425e15 uses different legacy inputs.",
)

# The theorem's MeV-line converts 49.53 meV as 49.53e-12 MeV. The correct
# conversion is 49.53e-9 MeV, or 49.53e-12 GeV.
sqrt_dm31_mev_correct = sqrt_dm31_ev * 1e-6
mr_correct_mev = (V_EW_GEV * 1000) ** 2 / sqrt_dm31_mev_correct
e_r_theorem_wrong = 65.189
v.check(
    "Theorem ER value after correcting meV-to-MeV conversion",
    PI + math.log(mr_correct_mev / M_E_MEV) / MU,
    e_r_theorem_wrong,
    abs_tol=0.05,
    detail="Expected fail: the theorem's 65.189 comes from a factor-1000 unit error.",
)
v.check("corrected E_R from listed inputs in later remark", e_r_from_listed_inputs, 56.485, abs_tol=0.01)
v.check("canonical P53/P60 E_R from table M_R", e_r_from_table_mr, 56.502, abs_tol=0.01)

gap_from_canonical = e_pl - e_r_canonical
gap_from_table_mr = e_pl - e_r_from_table_mr
gap_from_listed_inputs = e_pl - e_r_from_listed_inputs
v.check(
    "full-precision gap from CODATA Planck mass and table M_R",
    gap_from_table_mr,
    11.594,
    abs_tol=0.004,
    detail="Expected fail: exact recomputation gives about 11.5874.",
)
v.record(
    "gap using listed v_EW and Delta m31 inputs",
    abs(gap_from_listed_inputs - 11.594) < 0.004,
    f"{gap_from_listed_inputs:.6f}",
    "11.594 +/- 0.004",
    "Expected fail: the stated PDG inputs give about 11.6064.",
)
v.check("canonical subtraction 68.096 - 56.502", 68.096 - 56.502, 11.594, abs_tol=1e-12)

delta_canonical = H_E6 - 11.594
v.check("delta from canonical gap", delta_canonical, 0.406, abs_tol=1e-12)
v.check("canonical fractional deviation percent", 100 * delta_canonical / H_E6, 3.38, rel=2e-3)

gap_uncert = DM31_FRAC_UNCERT / (2 * MU)
v.check("gap uncertainty propagated from 1.4% Delta m31", gap_uncert, 0.0088, rel=5e-3)
v.check(
    "abstract gap uncertainty +/-0.004",
    gap_uncert,
    0.004,
    rel=1e-1,
    detail="Expected fail: the paper's own propagation gives about +/-0.0088.",
)
v.check("sigma exclusion using propagated uncertainty", delta_canonical / gap_uncert, 46.0, rel=2e-2)

v.check(
    "reduced Planck spectral coordinate",
    e_pl_reduced,
    66.476,
    abs_tol=0.01,
    detail="Expected fail: M_Pl/sqrt(8*pi) maps to about 66.060, not 66.476.",
)
v.check(
    "reduced Planck gap using canonical E_R",
    e_pl_reduced - e_r_canonical,
    9.974,
    abs_tol=0.02,
    detail="Expected fail: the corrected reduced-Planck gap is about 9.558.",
)

coupling_change = 38 * MU * H_E6 / (2 * PI)
v.check("E6 one-loop coupling change", coupling_change, 57.6, rel=1e-3)
v.check("coupling factor vs alpha_GUT^-1 about 25", coupling_change / 25, 2.3, rel=2e-2)
v.check("E6 level-1 WZW central charge", 78 / (1 + H_E6), 6.0, rel=1e-12)

e_r_conj = 68.096 - H_E6
mr_conj_from_e = mass_gev_from_e(e_r_conj)
mr_conj_from_pl = M_PL_GEV * math.exp(-MU * H_E6)
v.check("E_R if exact E6 gap uses paper canonical E_Pl", e_r_conj, 56.096, abs_tol=1e-12)
v.check("M_R from exact E6 gap via E_R", mr_conj_from_e, 8.95e14, rel=5e-3)
v.check("M_R from M_Pl*exp(-MU*12)", mr_conj_from_pl, 8.95e14, rel=5e-3)
mnu_conj_mev = V_EW_GEV**2 / mr_conj_from_pl * 1e12
v.check("m_nu3 from exact E6 gap", mnu_conj_mev, 67.7, rel=5e-3)
v.check("oscillation m_nu3 central value", sqrt_dm31_ev * 1000, 49.53, rel=1e-4)
v.check("exact-conjecture overestimate percent", 100 * (mnu_conj_mev / (sqrt_dm31_ev * 1000) - 1), 36.8, rel=3e-2)
v.check("exact-conjecture sigma exclusion", (mnu_conj_mev - 49.53) / 0.35, 52.0, rel=5e-2)

epsilon = 0.406
candidate_values = {
    "MU*pi/2": MU * PI / 2,
    "(1-MU)*pi/2": (1 - MU) * PI / 2,
    "MU*theta_C": MU * THETA_C,
    "pi/mu0": PI / MU0,
    "(pi^2-9)/MU": (PI**2 - 9) / MU,
    "1/pi^2": 1 / PI**2,
    "MU^2": MU**2,
    "(1-MU)^2/MU": (1 - MU) ** 2 / MU,
    "4*pi/mu0": 4 * PI / MU0,
    "sin^2(theta_C)": LAMBDA**2,
    "ln(1+MU)": math.log(1 + MU),
    "pi/(4*h_E6)": PI / (4 * H_E6),
    "(1-MU)*pi": (1 - MU) * PI,
}
claimed_candidates = {
    "MU*pi/2": 1.246,
    "(1-MU)*pi/2": 0.325,
    "MU*theta_C": 0.178,
    "pi/mu0": 0.0229,
    "(pi^2-9)/MU": 1.097,
    "1/pi^2": 0.1013,
    "MU^2": 0.6294,
    "(1-MU)^2/MU": 0.0539,
    "4*pi/mu0": 0.0917,
    "sin^2(theta_C)": 0.04952,
    "ln(1+MU)": 0.5848,
    "pi/(4*h_E6)": 0.0654,
    "(1-MU)*pi": 0.6493,
}
for name, claimed in claimed_candidates.items():
    rel = 2e-3
    if name == "ln(1+MU)":
        rel = 2e-3
    v.check(f"epsilon candidate {name}", candidate_values[name], claimed, rel=rel)

v.check(
    "fractional shortfall candidate MU^10",
    MU**10,
    0.0938,
    rel=1e-2,
    detail="Expected fail: 0.7933^10 is about 0.0988.",
)
v.record(
    "pi/mu0 is not alpha^-1",
    abs((PI / MU0) - MU0) < 1e-6,
    f"pi/mu0={PI / MU0:.6f}; alpha^-1=mu0={MU0:.6f}",
    "pi/mu0 = alpha^-1",
    "Expected fail: pi/mu0 is pi*alpha, not alpha inverse.",
)

closest = min(abs(value / epsilon - 1) for value in candidate_values.values())
v.record("no listed epsilon candidate within 10%", closest > 0.10, f"closest relative error {closest:.3f}", ">0.10")

sys.exit(v.summary())