#!/usr/bin/env python3
"""
verify_P076.py -- Addendum 76: OP-B-Mass CGF / spectral-weight closure.

This verifier checks the CGF algebra and numerical appendix values in
76_Addendum_OPBMassClosure.tex.  The basic rho-derived identities are solid:
CGF(0)=ln(mu0), CGF'(0)=MU, the log-charge gap is ln(mu0)/MU, and the mass
map sends that exact gap to mu0.  The flagged issues are stale precise table
values, a slightly off entropy computation, and proof/status claims around
lambda=1 uniqueness and the polynomial/transcendence bridge.
"""

from __future__ import annotations

import math
import sys
from pathlib import Path

import mpmath as mp

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


mp.mp.dps = 60
v = Verifier("P076 -- OP-B-Mass CGF Closure")

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

PI = mp.pi
MU0 = 4 * PI**3 + PI**2 + PI
MU1 = 16 * PI**3 / 5 + 3 * PI**2 / 4 + 2 * PI / 3
MU = MU1 / MU0
MU2 = 2 * PI / 4 + 3 * PI**2 / 5 + 16 * PI**3 / 6
CGF0 = mp.log(MU0)
CGF1 = MU
KAPPA2 = MU2 / MU0 - MU**2
GAP = CGF0 / CGF1
E_EDGE = PI
E_BDY = PI**2
E_BULK = 4 * PI**3
E_CHARM = E_EDGE + E_BDY
E_TOP = E_CHARM + GAP


def f(x: mp.mpf) -> float:
    return float(x)


def rho(x: mp.mpf) -> mp.mpf:
    return 2 * PI * x + 3 * PI**2 * x**2 + 16 * PI**3 * x**3


def mgf(t: mp.mpf) -> mp.mpf:
    return mp.quad(lambda x: mp.e ** (t * x) * rho(x), [0, 1])


def mgf_formula(t: mp.mpf) -> mp.mpf:
    return (
        2 * PI / t**2 * (mp.e**t * (t - 1) + 1)
        + 3 * PI**2 / t**3 * (mp.e**t * (t**2 - 2 * t + 2) - 2)
        + 16 * PI**3 / t**4 * (mp.e**t * (t**3 - 3 * t**2 + 6 * t - 6) + 6)
    )


v.check("mu0", f(MU0), 137.036304, rel=5e-9)
v.check("mu1", f(MU1), 108.716684, rel=5e-9)
v.check("MU", f(MU), 0.793342, rel=5e-7)
v.check("E_edge", f(E_EDGE), math.pi, rel=1e-12)
v.check("E_boundary", f(E_BDY), math.pi**2, rel=1e-12)
v.check("E_bulk", f(E_BULK), 124.025, rel=1e-5)
v.check("E_c exact", f(E_CHARM), 13.011378, rel=5e-6, detail="Expected fail: exact pi+pi^2 is about 13.011197.")

v.check("MGF(0)", f(MU0), 137.036304, rel=5e-9)
v.check("CGF(0) rounded in prose", f(CGF0), 4.920, rel=2e-4)
v.check(
    "CGF(0) appendix value",
    f(CGF0),
    4.919716,
    rel=2e-5,
    detail="Expected fail: the precise appendix value is stale; ln(mu0)≈4.920246.",
)
v.check("CGF'(0)", f(CGF1), 0.793342, rel=5e-7)
v.check("mu2", f(MU2), 90.17, rel=1e-4)
v.check("mu2/mu0", f(MU2 / MU0), 0.658, rel=8e-4)
v.check("MU^2", f(MU**2), 0.629, rel=8e-4)
v.check("kappa2", f(KAPPA2), 0.029, rel=2e-2)

v.check("MGF explicit formula at t=0.3", f(mgf_formula(mp.mpf("0.3"))), f(mgf(mp.mpf("0.3"))), rel=1e-12)
v.check("CGF ratio rounded in prose", f(GAP), 6.201, rel=2e-4)
v.check(
    "CGF ratio appendix value",
    f(GAP),
    6.201027,
    rel=5e-5,
    detail="Expected fail: exact ln(mu0)/MU≈6.201921.",
)
v.check("log-charge formula", f(MU0 * CGF0 / MU1), f(GAP), rel=1e-12)
v.check("saturation equation exp(MU gap)", f(mp.e ** (MU * GAP)), f(MU0), rel=1e-12)
v.check("predicted top energy exact", f(E_TOP), 19.212405, rel=2e-5, detail="Expected fail: exact Ec+ln(mu0)/MU≈19.213118.")

printed_gap = mp.mpf("19.212405") - mp.mpf("13.011378")
printed_ratio = mp.e ** (MU * printed_gap)
v.check(
    "printed appendix E_c/E_t reproduce top-charm ratio",
    f(printed_ratio),
    137.036304,
    rel=5e-4,
    detail="Expected fail: using the printed E_c and E_t gives about 136.939, not mu0.",
)
v.check("PDG ratio residual using 137.06", 100 * f((MU0 - mp.mpf("137.06")) / mp.mpf("137.06")), -0.017, rel=3e-2)

v.check("bulk spectral fraction", f(E_BULK / MU0), 0.9051, rel=7e-5)
v.check("charm spectral fraction", f(E_CHARM / MU0), 0.0949, rel=5e-4)
v.check("bulk-only gap", f(mp.log(E_BULK) / MU), 6.077, rel=2e-4)
v.check("bulk-only top energy", f(E_CHARM + mp.log(E_BULK) / MU), 19.089, rel=1e-4)
v.check("bulk-only ratio", f(mp.e ** (MU * (mp.log(E_BULK) / MU))), 124.0, rel=3e-4)
v.check("electron-baseline top energy contradiction", f(E_EDGE + GAP), 9.344, rel=1e-4)
v.check("cubic norm", f(E_EDGE * E_BDY * E_BULK), f(4 * PI**6), rel=1e-12)

entropy_integral = mp.quad(lambda x: rho(x) * mp.log(rho(x)), [0, 1])
entropy = CGF0 - entropy_integral / MU0
v.check("entropy integral exact numeric", f(entropy_integral), 756.0, rel=5e-4, detail="Expected fail: direct quadrature gives about 756.819.")
v.check(
    "normalized differential entropy",
    f(entropy),
    -0.597,
    rel=5e-3,
    detail="Expected fail: using the exact integral gives about -0.6025.",
)
v.check("KL divergence from uniform", f(-entropy), 0.597, rel=5e-3, detail="Expected fail: direct quadrature gives about 0.6025.")

z_values = [mp.quad(lambda x, lam=lam: mp.e ** (lam * 0 * x) * rho(x), [0, 1]) for lam in [0, 1, 2]]
v.record(
    "lambda=1 is selected by Z(lambda,0)",
    not all(abs(z - MU0) < mp.mpf("1e-45") for z in z_values),
    computed=", ".join(f"lambda {lam}: {float(z):.6f}" for lam, z in zip([0, 1, 2], z_values)),
    claimed="lambda=1 saturation is distinguished at t=0",
    detail="Expected proof-audit fail: at t=0 the exponential is 1, so Z(lambda,0)=mu0 for every lambda.",
)
v.record(
    "lambda=1 saturation is unique in lambda",
    False,
    computed="the equation e^{MU Delta}=Z(lambda,0)=mu0 has the same Delta for every lambda",
    claimed="lambda=1 condition is the unique saturation",
    detail="Expected proof-audit fail: uniqueness holds for Delta after choosing mu0, not for lambda.",
)
v.record(
    "ln-bridge transcendence proof is rigorous as written",
    False,
    computed="proof asserts ln(Q(pi)) is not in Q(pi) and cites transcendence of e",
    claimed="no polynomial bridge established",
    detail="Expected proof-audit fail: the stated transcendence argument does not prove the theorem for Q(pi).",
)
v.record(
    "SWS residual is explicitly still open",
    "Top as $\\lambda=1$ saturation point (SWS)" in TEX and "\\textbf{Open}" in TEX,
    computed="status table marks SWS open",
    claimed="conditional closure only",
)

sys.exit(v.summary())