"""
verify_P121.py
Standalone verification of Addendum P121:
  "Strong Coupling αs at NLO: b₁ = 116/3 from the J₃(𝕆) Two-Loop Peirce-Trace"

Paper claims:
  - b₀ = 23/3  (one-loop β-function coefficient, nf=5, SU(3))
  - b₁ = 116/3 ≈ 38.67  (two-loop coefficient from J₃(𝕆) Peirce Casimirs)
  - b₁/b₀ = 116/23 ≈ 5.043
  - mconf = π · ALPHA_INV · m_e ≈ 219.99 MeV  (BREATH_PERIOD confinement scale)
  - ln(MZ/mconf) ≈ 6.027
  - δ ≈ 0.03402  (fractional NLO shift)
  - δαs ≈ −0.0040  (absolute NLO shift)
  - αs_NLO(MZ) = 0.1186 · (1 − δ) ≈ 0.1146  (PRIMARY CLAIM)
  - Bracketing: αs_NLO = 0.1146 < αs_PDG = 0.1179 < αs_LO = 0.1186

NLO formula (paper eq. 13):
  αs_NLO = αs_LO · [1 − (b₁/b₀) · αs_LO² · ln(MZ/mconf) / (4π)]

β-function coefficients (paper eqs. 2–3, SU(Nc) with nf massless quarks):
  b₀ = (11/3)·CA − (4/3)·TF·nf
  b₁ = (34/3)·CA² − ((20/3)·CA + 4·CF)·TF·nf

Casimirs (Proposition 2.1, J₃(𝕆) Peirce decomposition):
  CA = 3,  CF = 4/3,  TF = 1/2

Copyright: Léon Fernando Vlegels. License: MIT.
"""

import math
import sys

try:
    from scipy.integrate import solve_ivp
    _SCIPY = True
except ImportError:
    _SCIPY = False

# ──────────────────────────────────────────────────────────────────
# 0. Constants
# ──────────────────────────────────────────────────────────────────

# TOE fine-structure inverse (kernel constant, P01/P36)
ALPHA_INV = 4 * math.pi**3 + math.pi**2 + math.pi   # ≈ 137.036
ALPHA     = 1.0 / ALPHA_INV
BREATH_PERIOD = math.pi * ALPHA_INV                  # ≈ 432 system units

# Physical constants
M_E_MEV  = 0.51099895   # electron mass (MeV)
MZ_MEV   = 91187.6      # Z boson mass (MeV)
AS_PDG   = 0.1179       # PDG αs(MZ) central value
AS_PDG_ERR = 0.0009     # PDG uncertainty (1σ)

# LO boundary condition from Addendum 106
AS_LO    = 0.1186       # multi-threshold 1-loop result
AS_BC    = math.sqrt(3) # αs(mconf) = cot(π/6) = √3  (G₂ root-length ratio)

# SU(3)_c Casimir invariants (Proposition 2.1)
CA  = 3        # adjoint Casimir of SU(3)
CF  = 4.0/3.0  # fundamental Casimir = (N²-1)/(2N)|N=3
TF  = 0.5      # trace normalisation = 1/2 for SU(N)
NF  = 5        # active flavours between mconf and MZ

# ──────────────────────────────────────────────────────────────────
# 1. β-function coefficient helpers
# ──────────────────────────────────────────────────────────────────

def compute_b0(ca=CA, tf=TF, nf=NF):
    """One-loop: b₀ = (11/3)·CA − (4/3)·TF·nf"""
    return (11.0/3.0)*ca - (4.0/3.0)*tf*nf

def compute_b1(ca=CA, cf=CF, tf=TF, nf=NF):
    """Two-loop: b₁ = (34/3)·CA² − ((20/3)·CA + 4·CF)·TF·nf"""
    gluon = (34.0/3.0) * ca**2
    quark = ((20.0/3.0)*ca + 4.0*cf) * tf * nf
    return gluon - quark

def nlo_shift(as_lo, b0, b1, ln_ratio):
    """
    Paper eq. 13: fractional NLO correction δ, so αs_NLO = αs_LO·(1−δ).
    δ = (b₁/b₀) · αs_LO² · ln(MZ/mconf) / (4π)
    """
    return (b1/b0) * as_lo**2 * ln_ratio / (4.0 * math.pi)

# ──────────────────────────────────────────────────────────────────
# 2. Full two-loop RGE (numerical cross-check only)
# ──────────────────────────────────────────────────────────────────

def rge_2loop(as_start, mu_start_mev, mu_end_mev, b0, b1):
    """
    Integrate the 2-loop RGE dαs/d(ln μ) = −(b₀/2π)·αs² − (b₁/4π²)·αs³
    from mu_start to mu_end using scipy solve_ivp.

    Note: This uses a SINGLE threshold (nf=5) throughout and starts from
    the TOE non-perturbative boundary αs(mconf) = √3, giving a different
    LO value (≈0.1261) than the multi-threshold A106 result (0.1186).
    The perturbative NLO formula in the paper is applied on top of the
    multi-threshold LO, so this ODE result is a consistency check on the
    sign and order-of-magnitude of the 2-loop correction only.
    """
    if not _SCIPY:
        return None

    t0 = math.log(mu_start_mev)
    t1 = math.log(mu_end_mev)

    def beta(t, y):
        a = y[0]
        return [-(b0/(2.0*math.pi))*a**2 - (b1/(4.0*math.pi**2))*a**3]

    sol = solve_ivp(beta, [t0, t1], [as_start], method='DOP853',
                    rtol=1e-10, atol=1e-12, dense_output=False)
    if sol.success:
        return sol.y[0, -1]
    return None

# ──────────────────────────────────────────────────────────────────
# 3. Main verification
# ──────────────────────────────────────────────────────────────────

PASS = FAIL = 0
_N = 0


def _mark(desc, cond):
    """Modern check line: increments PASS/FAIL counters."""
    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


def run():
    failures = []

    def check(name, got, expected, tol, *, fmt=".6g"):
        ok = abs(got - expected) <= tol
        _mark(name, ok)
        print(f"         got={got:{fmt}}  expected={expected:{fmt}}  tol=±{tol:.2g}")
        if not ok:
            failures.append((name, got, expected, tol))

    # ──────────────────────────────────────────────
    print("=" * 65)
    print("P121 Strong Coupling αs at NLO — Verification")
    print("=" * 65)
    print()

    # ── A. TOE constants ──────────────────────────
    print("── A. TOE constants ────────────────────────────────────────")
    print(f"  ALPHA_INV  = {ALPHA_INV:.8f}  (paper: ~137.036)")
    print(f"  ALPHA      = {ALPHA:.10f}")
    print(f"  BREATH_PERIOD = π · ALPHA_INV = {BREATH_PERIOD:.4f}  (≈432)")
    print()

    # ── B. Confinement scale ──────────────────────
    print("── B. Confinement scale ────────────────────────────────────")
    m_conf = math.pi * ALPHA_INV * M_E_MEV   # = BREATH_PERIOD · m_e
    print(f"  mconf = π · ALPHA_INV · m_e = {m_conf:.4f} MeV")
    print(f"  (paper claims ≈ 219.99 MeV;  BREATH_PERIOD = {BREATH_PERIOD:.4f})")
    check("mconf [MeV]", m_conf, 219.99, 0.05)
    print()

    # ── C. Casimir invariants ─────────────────────
    print("── C. Casimir invariants (SU(3)_c ⊂ G₂ ⊂ F₄) ────────────")
    # CA = N_c = 3  (adjoint quadratic Casimir of SU(3))
    CA_check = 3
    # CF = (N²-1)/(2N) for N=3
    CF_check = (3**2 - 1) / (2.0 * 3)
    # TF = 1/2 (standard generator normalisation)
    TF_check = 0.5
    print(f"  CA = {CA_check}  (paper: 3)")
    print(f"  CF = (N²−1)/(2N)|N=3 = {CF_check:.4f}  (paper: 4/3 = {4/3:.4f})")
    print(f"  TF = {TF_check}  (paper: 1/2)")
    check("CA",    CA_check,  3,       1e-12)
    check("CF",    CF_check,  4.0/3.0, 1e-12)
    check("TF",    TF_check,  0.5,     1e-12)
    print()

    # ── D. β-function coefficients ────────────────
    print("── D. β-function coefficients ──────────────────────────────")
    b0 = compute_b0()
    b0_expected = 23.0/3.0
    b1_gluon = (34.0/3.0) * CA**2
    b1_quark = -((20.0/3.0)*CA + 4.0*CF) * TF * NF
    b1 = b1_gluon + b1_quark
    b1_expected = 116.0/3.0
    ratio = b1 / b0
    ratio_expected = 116.0/23.0

    print(f"  b₀ = (11/3)·CA − (4/3)·TF·nf  = {b0:.8f}  (expected 23/3 = {b0_expected:.8f})")
    print(f"  b₁_gluon = (34/3)·CA²          = {b1_gluon:.8f}  (expected 102)")
    print(f"  b₁_quark = −[(20/3)·CA+4·CF]·TF·nf = {b1_quark:.8f}  (expected −190/3 = {-190/3:.8f})")
    print(f"  b₁ = b₁_gluon + b₁_quark       = {b1:.8f}  (expected 116/3 = {b1_expected:.8f})")
    print(f"  b₁/b₀                           = {ratio:.8f}  (expected 116/23 = {ratio_expected:.8f})")
    check("b₀",    b0,    b0_expected,    1e-12)
    check("b₁",    b1,    b1_expected,    1e-12)
    check("b₁/b₀", ratio, ratio_expected, 1e-12)

    # Cross-check: user-convention coefficients.
    # Alt convention writes dαs/d(lnμ) = −2b̃₀αs² − 4b̃₁αs³
    # so the paper's coefficients satisfy b₀_paper = 4π·b̃₀,  b₁_paper = 16π²·b̃₁
    # where b̃₀ = (33−2nf)/(12π)  and  b̃₁ = (153−19nf)/(24π²).
    #
    # Derivation: paper eq.1 uses dαs/dt = −(b₀/2π)αs² − (b₁/4π²)αs³  (t=lnμ).
    # Alt convention: same equation written as −2b̃₀αs² − 4b̃₁αs³
    #   ⟹ b₀/(2π) = 2b̃₀  ⟹  b₀ = 4π·b̃₀
    #   ⟹ b₁/(4π²) = 4b̃₁  ⟹  b₁ = 16π²·b̃₁
    beta0_alt = (33.0 - 2.0*NF) / (12.0 * math.pi)       # b̃₀
    beta1_alt = (153.0 - 19.0*NF) / (24.0 * math.pi**2)  # b̃₁
    b0_from_alt = 4.0 * math.pi * beta0_alt               # = b₀_paper
    b1_from_alt = 16.0 * math.pi**2 * beta1_alt           # = b₁_paper
    print(f"  Cross-check (alt convention dαs/dt = −2b̃₀αs²−4b̃₁αs³, nf=5):")
    print(f"    b̃₀ = (33−2nf)/(12π)          = {beta0_alt:.8f}")
    print(f"    b̃₁ = (153−19nf)/(24π²)       = {beta1_alt:.8f}")
    print(f"    4π·b̃₀  → b₀_paper            = {b0_from_alt:.8f}")
    print(f"    16π²·b̃₁ → b₁_paper           = {b1_from_alt:.8f}")
    check("b₀ from alt convention",  b0_from_alt, b0_expected, 1e-10)
    check("b₁ from alt convention",  b1_from_alt, b1_expected, 1e-10)
    print()

    # ── E. NLO shift ──────────────────────────────
    print("── E. NLO shift (perturbative, paper eq. 13) ───────────────")
    ln_ratio = math.log(MZ_MEV / m_conf)
    delta    = nlo_shift(AS_LO, b0, b1, ln_ratio)
    delta_as = -AS_LO * delta
    as_nlo   = AS_LO * (1.0 - delta)

    print(f"  αs_LO (Addendum 106)            = {AS_LO:.4f}")
    print(f"  MZ/mconf                         = {MZ_MEV/m_conf:.3f}")
    print(f"  ln(MZ/mconf)                     = {ln_ratio:.6f}  (paper: 6.027)")
    print(f"  δ = (b₁/b₀)·αs_LO²·ln(MZ/mc)/(4π) = {delta:.6f}  (paper: 0.03402)")
    print(f"  δαs = −αs_LO · δ                = {delta_as:.6f}  (paper: −0.004)")
    print(f"  αs_NLO(MZ) = αs_LO·(1−δ)       = {as_nlo:.4f}  (paper: 0.1146)")

    check("ln(MZ/mconf)",   ln_ratio, 6.027,  0.001)
    check("δ (frac shift)", delta,    0.03402, 0.0001)
    check("δαs",            delta_as, -0.004,  0.0001)
    check("αs_NLO(MZ)",     as_nlo,   0.1146,  0.0001, fmt=".6f")  # PRIMARY CLAIM
    print()

    # ── F. Bracketing ─────────────────────────────
    print("── F. Bracketing: αs_NLO < αs_PDG < αs_LO ────────────────")
    bracket_ok = as_nlo < AS_PDG < AS_LO
    sigma_nlo  = (as_nlo - AS_PDG) / AS_PDG_ERR
    sigma_lo   = (AS_LO  - AS_PDG) / AS_PDG_ERR
    print(f"  αs_NLO = {as_nlo:.4f} < αs_PDG = {AS_PDG:.4f} < αs_LO = {AS_LO:.4f}  → {bracket_ok}")
    print(f"  NLO significance: {sigma_nlo:.2f}σ  (paper: −3.68σ)")
    print(f"  LO  significance: {sigma_lo:+.2f}σ  (paper: +0.80σ)")
    _mark("PDG bracketed between NLO and LO", bracket_ok)
    if not bracket_ok:
        failures.append(("bracketing", as_nlo, AS_PDG, 0))
    check("NLO sigma [σ]", sigma_nlo, -3.68, 0.05)
    check("LO  sigma [σ]", sigma_lo,  +0.80, 0.05)
    print()

    # ── G. Full 2-loop RGE numerical integration ──
    print("── G. Full 2-loop RGE numerical integration (informational) ")
    if _SCIPY:
        # Single-threshold nf=5 from boundary αs(mconf)=√3.
        # This gives a different LO (≈0.1261) from A106 multi-threshold (0.1186),
        # so the 2-loop result here is a cross-check on sign/direction only.
        as_2loop_single = rge_2loop(AS_BC, m_conf, MZ_MEV, b0, b1)
        as_1loop_single = rge_2loop(AS_BC, m_conf, MZ_MEV, b0, 0.0)
        if as_2loop_single is not None:
            print(f"  Boundary: αs(mconf={m_conf:.2f} MeV) = √3 = {AS_BC:.6f}")
            print(f"  1-loop ODE result (nf=5 single-threshold): αs(MZ) = {as_1loop_single:.4f}")
            print(f"  2-loop ODE result (nf=5 single-threshold): αs(MZ) = {as_2loop_single:.4f}")
            print(f"  2-loop shift vs 1-loop: Δαs = {as_2loop_single - as_1loop_single:+.4f}")
            nlo_sign_ok = as_2loop_single < as_1loop_single
            print(f"  NLO moves αs downward (correct sign): {nlo_sign_ok}")
            _mark("2-loop shift direction consistent with paper", nlo_sign_ok)
            if not nlo_sign_ok:
                failures.append(("2-loop RGE sign", as_2loop_single, as_1loop_single, 0))
        else:
            print("  [SKIP] scipy solve_ivp failed; skipping numerical RGE")
    else:
        print("  [SKIP] scipy not available; install with: pip install scipy")
    print()

    # ──────────────────────────────────────────────
    # Summary
    # ──────────────────────────────────────────────
    if failures:
        for name, got, expected, tol in failures:
            print(f"  ✗ {name}: got {got:.6g}, expected {expected:.6g}, tol ±{tol:.2g}")
    else:
        print("Primary claim: αs_NLO(MZ) = "
              f"{as_nlo:.4f}  (paper: 0.1146)  ✓")
        print(f"  b₁ = 116/3 = {116/3:.4f}  ✓")
        print(f"  δαs = {delta_as:.4f}  (paper: −0.0040)  ✓")
        print(f"  PDG bracketed: {as_nlo:.4f} < {AS_PDG:.4f} < {AS_LO:.4f}  ✓")

    print(f"\n{'='*60}\nRESULT: {PASS} PASS / {FAIL} FAIL")
    sys.exit(0 if FAIL == 0 else 1)


if __name__ == "__main__":
    run()