#!/usr/bin/env python3
"""
verify_P098.py -- Addendum 98: Bryant-Fano NLO loop and Jarlskog protection.

This verifier checks the two-channel cancellation arithmetic, the updated
Wolfenstein/Jarlskog table, and selected proof-status claims in
98_Addendum_BryantFanoNLO.tex.

The high-level difference-of-squares arithmetic is reproducible if one assumes
opposite +/-4 lambda^2/5 loop factors. The flagged issues are that the paper's
own Moufang-loop correction formula first gives only +lambda^2/N, so N=5 would
not cancel the G2 loop; the proof obtains -1/4 on the CP direction before
reinterpreting it as a positive correction; the abstract uses a sqrt-style
1-8 lambda^4/25 factor for A while the body/table use 1-16 lambda^4/25; and
the exact Phi/J route remains stale relative to the printed trigonometric formula.
"""

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("P098 -- Bryant-Fano NLO Jarlskog")

PI = math.pi
LAM = math.sin(PI / 14)
LAM2 = LAM**2
LAM4 = LAM2**2
LAM6 = LAM**6
A_LO = math.cos(PI / 14) ** 2 * math.cos(2 * PI / 14)
PHI_EXACT = math.cos(PI / 14) * math.cos(2 * PI / 14) * math.cos(5 * PI / 28)
DELTA = 31 * PI / 84
SIN_DELTA = math.sin(DELTA)
COS_DELTA = math.cos(DELTA)
ETA_EXACT = 0.5 * PHI_EXACT * SIN_DELTA
RHO_EXACT = 0.5 * PHI_EXACT * COS_DELTA
J_LO_EXACT = A_LO**2 * LAM6 * ETA_EXACT

J_LO_ROUNDED = 3.064e-5
J_PDG = 3.08e-5
J_SIGMA = 0.15e-5
PHI_STALE = 0.7429
ETA_STALE = 0.3439
RHO_STALE = 0.1484
A_STALE = 0.8563


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


v.check("lambda = sin(pi/14)", LAM, 0.222521, rel=5e-6)
v.check("lambda^2", LAM2, 0.04951, rel=2e-4)
v.check("lambda^4", LAM4, 0.002451, rel=5e-4)
v.check("lambda^6", LAM6, 1.2148e-4, rel=8e-4)
v.check("A_LO", A_LO, 0.8563, rel=1e-4)
v.check(
    "Phi_LO from exact formula",
    PHI_EXACT,
    0.7429,
    rel=5e-4,
    detail="Expected fail: the printed exact formula gives about 0.743745.",
)
v.check("sin(delta_0)", SIN_DELTA, 0.9162, rel=5e-4)
v.check("cos(delta_0)", COS_DELTA, 0.4008, rel=3e-3)
v.check(
    "J_LO from exact stated formulas",
    J_LO_EXACT,
    3.064e-5,
    rel=5e-4,
    detail="Expected fail: exact Phi-derived eta gives about 3.035e-5; 3.064e-5 is the stale rounded path.",
)

g2_factor = 1 - 4 * LAM2 / 5
fano_factor_assumed = 1 + 4 * LAM2 / 5
combined_block = g2_factor * fano_factor_assumed
j_ratio_exact = combined_block**3
j_ratio_first_order = 1 - 48 * LAM4 / 25
v.check("G2 loop factor 1 - 4 lambda^2/5", g2_factor, 0.96039, rel=5e-6)
v.check("assumed Fano loop factor 1 + 4 lambda^2/5", fano_factor_assumed, 1.03961, rel=5e-6)
v.check("combined block factor", combined_block, 1 - 16 * LAM4 / 25, rel=1e-12)
v.check("16 lambda^4 / 25", 16 * LAM4 / 25, 0.0015685, rel=5e-4)
v.check("48 lambda^4 / 25", 48 * LAM4 / 25, 0.0047059, rel=5e-4)
v.check("J ratio exact (1 - 16 lambda^4/25)^3", j_ratio_exact, 0.99529, rel=5e-5)
v.check("J ratio first-order lambda^4", j_ratio_first_order, 0.995294, rel=5e-5)
v.check("rounded-path J_NLO two-channel", J_LO_ROUNDED * j_ratio_exact, 3.050e-5, rel=2e-4)
v.check("rounded-path J residual percent", residual_pct(J_LO_ROUNDED * j_ratio_exact, J_PDG), -0.97, rel=3e-2)
v.check("rounded-path J sigma", abs(J_LO_ROUNDED * j_ratio_exact - J_PDG) / J_SIGMA, 0.20, rel=2e-2)

v.check(
    "exact-formula J_NLO two-channel vs claimed",
    J_LO_EXACT * j_ratio_exact,
    3.050e-5,
    rel=5e-3,
    detail="Expected fail: exact Phi-derived J_LO gives about 3.020e-5 after the same ratio.",
)
v.check(
    "exact-formula J residual percent",
    residual_pct(J_LO_EXACT * j_ratio_exact, J_PDG),
    -0.97,
    rel=5e-2,
    detail="Expected fail: exact formula path gives about -1.94%, not -0.97%.",
)

# Internal formula consistency for the Fano correction.
L_FANO = 1 / 4
SIGMA = 4 * LAM2
N_FANO = 5
fano_delta_from_displayed_formula = L_FANO * SIGMA / N_FANO
fano_delta_needed = 4 * LAM2 / 5
v.check(
    "displayed Fano correction formula gives +4 lambda^2/5 with N=5",
    fano_delta_from_displayed_formula,
    fano_delta_needed,
    rel=1e-6,
    detail="Expected fail: (1/4)*(4 lambda^2)/5 = lambda^2/5, a factor of four too small.",
)
v.check(
    "N_Fano required by displayed +lambda^2/N formula",
    LAM2 / fano_delta_needed,
    5,
    rel=1e-6,
    detail="Expected fail: the displayed +lambda^2/N formula would require N=5/4 to cancel -4 lambda^2/5.",
)
v.record(
    "Moufang eigenvalue is +1/4 on the CP e7 direction",
    False,
    "The proof computes M_a(E23(e7)) = -1/4 E23(e7) before reinterpreting the correction sign.",
    "+1/4 eigenvalue",
    "Expected proof-audit fail: the positive correction is an extra sign convention, not the eigenvalue just computed on the CP direction.",
)
v.record(
    "Moufang and Peirce maps are simultaneously equal and independent",
    False,
    "Definition states a o (X o a) = (a o X) o a, while the two-channel proof says these maps are distinct in J3(O).",
    "consistent map identities",
    "Expected proof-audit fail: the equality in the definition conflicts with the independence proof.",
)
v.record(
    "N_Fano=5 is derived from first principles",
    False,
    "The paper derives N=5 from the cancellation condition and marks the root-system derivation as OP-P98-1.",
    "first-principles derivation",
    "Expected proof-audit fail: the value is conditional/self-consistency based, not first-principles proved.",
)

# Wolfenstein table checks.
body_block = 1 - 16 * LAM4 / 25
sqrt_style = 1 - 8 * LAM4 / 25
v.check(
    "abstract A_NLO factor 1 - 8 lambda^4/25 matches body/table factor",
    sqrt_style,
    body_block,
    rel=1e-6,
    detail="Expected fail: the abstract uses the square-root expansion; the body/table use 1 - 16 lambda^4/25.",
)
v.check("two-channel A table", A_STALE * body_block, 0.8550, rel=1e-4)
v.check("two-channel Phi table", PHI_STALE * body_block, 0.7417, rel=5e-5)
v.check("two-channel eta table", ETA_STALE * body_block, 0.3434, rel=2e-4)
v.check("two-channel rho table", RHO_STALE * body_block, 0.1482, rel=3e-4)
v.check("G2-only rounded J table", J_LO_ROUNDED * g2_factor**3, 2.714e-5, rel=5e-4)
v.check(
    "G2-only exact-formula J table",
    J_LO_EXACT * g2_factor**3,
    2.714e-5,
    rel=5e-3,
    detail="Expected fail: exact formula path gives about 2.688e-5.",
)

eta_lo_resid = abs(residual_pct(ETA_STALE, 0.3441))
eta_two_resid = abs(residual_pct(ETA_STALE * body_block, 0.3441))
v.record(
    "two-channel NLO parameters are closer to PDG for eta",
    eta_two_resid < eta_lo_resid,
    f"LO eta residual={eta_lo_resid:.4f}%, two-channel={eta_two_resid:.4f}%",
    "two-channel closer",
    "Expected fail: eta moves slightly farther from the PDG central value.",
)
v.record(
    "16 lambda^4/25 is sub-permil",
    16 * LAM4 / 25 < 0.001,
    f"{1000 * 16 * LAM4 / 25:.3f} per mille",
    "<1 per mille",
    "Expected wording fail: the correction is about 1.57 per mille, i.e. sub-percent but not sub-permil.",
)

sys.exit(v.summary())
