"""
verify_P132.py — Numerical verification of Addendum 132
Structural Origin of the Matching Scale μ_match ≈ 706 MeV

Checks every numerical claim in the paper.
Run with:  python Lumen/corpus/addenda/verify/verify_P132.py
"""

import math
import sys

# ── TOE constants ────────────────────────────────────────────────────────────
ALPHA_INV = 4 * math.pi**3 + math.pi**2 + math.pi   # ≈ 137.036
ALPHA     = 1 / ALPHA_INV
BREATH_PERIOD = math.pi * ALPHA_INV                   # ≈ 432

# ── Physical inputs ──────────────────────────────────────────────────────────
M_E_MEV   = 0.511          # electron mass, MeV
N_F       = 5              # active quark flavours (P127/P132 convention)
C_A       = 3              # QCD colour factor
T_F       = 0.5            # quark Dynkin index

# ── G₂ root geometry (all exact / algebraic) ─────────────────────────────────
RHO_G2    = math.sqrt(3)   # G₂ long-to-short root-length ratio
PHI_POS   = 6              # |Φ⁺(G₂)|, positive-root count

# ── J₃(𝕆) / one-loop QCD coefficient (n_f = 5) ──────────────────────────────
# b₀ = (11 C_A − 4 T_F n_f) / 3  [standard MSbar one-loop β₀]
B0        = (11 * C_A - 4 * T_F * N_F) / 3   # = 23/3

# ── Tolerance helpers ────────────────────────────────────────────────────────
PASS = FAIL = 0
_N = 0

def _mark(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

def _next():
    global _N
    _N += 1
    return _N

def check(label, got, expected, tol=1e-3, unit=""):
    err = abs(got - expected) / abs(expected)
    ok  = err < tol
    return _mark(_next(),
                 f"{label}: got {got:.6g}{unit}  expected ≈ {expected:.6g}{unit}  (rel err {err:.2e})",
                 ok)


def check_exact(label, got, expected, tol=1e-12):
    """Exact algebraic identity — tight tolerance."""
    err = abs(got - expected)
    ok  = err < tol
    return _mark(_next(),
                 f"{label}: got {got!r}  expected {expected!r}  (abs err {err:.2e})",
                 ok)


# ════════════════════════════════════════════════════════════════════════════
print("P132 — Structural Origin of the Matching Scale μ_match ≈ 706 MeV")

# ────────────────────────────────────────────────────────────────────────────
print("S1  TOE constants")
# ────────────────────────────────────────────────────────────────────────────
check_exact("α⁻¹ = 4π³+π²+π", ALPHA_INV, 4*math.pi**3 + math.pi**2 + math.pi)
check("α⁻¹ (numeric)", ALPHA_INV, 137.036, tol=1e-3)

# ────────────────────────────────────────────────────────────────────────────
print("S2  G₂ root geometry identities")
# ────────────────────────────────────────────────────────────────────────────
# 2 ρ_G₂² = 6 = |Φ⁺(G₂)|
check_exact("2·ρ_G₂² = |Φ⁺(G₂)| = 6", 2 * RHO_G2**2, float(PHI_POS))

# ρ_G₂(2ρ_G₂ − 1) = 6 − √3  (algebraic simplification, §2 eq.(5))
lhs = RHO_G2 * (2*RHO_G2 - 1)
rhs = float(PHI_POS) - RHO_G2
check_exact("ρ_G₂(2ρ_G₂−1) = 6−√3", lhs, rhs)

# α_s at confinement scale: cot(π/6) = √3
alpha_s_conf = 1 / math.tan(math.pi / 6)
check_exact("α_s(m_conf) = cot(π/6) = √3", alpha_s_conf, math.sqrt(3))

# α_s at matching scale: sin(π/6) = 1/2
alpha_s_match = math.sin(math.pi / 6)
check_exact("α_s(μ_match) = sin(π/6) = 1/2", alpha_s_match, 0.5)

# ────────────────────────────────────────────────────────────────────────────
print("S3  J₃(𝕆) one-loop coefficient b₀ (n_f = 5)")
# ────────────────────────────────────────────────────────────────────────────
# Gluonic: 11 C_A / 3 = 11
gluon_part = 11 * C_A / 3
check_exact("gluonic contribution 11·C_A/3 = 11", gluon_part, 11.0)

# Quark-loop subtraction: −4 T_F n_f / 3 = −10/3
quark_part = -4 * T_F * N_F / 3
check_exact("quark subtraction −4·T_F·n_f/3 = −10/3", quark_part, -10.0/3)

# b₀ = 23/3
check_exact("b₀ = 23/3", B0, 23.0/3)
print(f"          [INFO] b₀ = {B0:.6f},  3·b₀ = {3*B0:.1f}")

# ────────────────────────────────────────────────────────────────────────────
print("S4  Confinement scale m_conf = π · α⁻¹ · m_e")
# ────────────────────────────────────────────────────────────────────────────
M_CONF_MEV = math.pi * ALPHA_INV * M_E_MEV
check("m_conf (MeV)", M_CONF_MEV, 219.991, tol=1e-3, unit=" MeV")
print(f"          [INFO] m_conf = {M_CONF_MEV:.4f} MeV")

# ────────────────────────────────────────────────────────────────────────────
print("S5  Exponent and scale factor")
# ────────────────────────────────────────────────────────────────────────────
# Correct exponent: 2π(6−√3)/23
exponent_correct = 2 * math.pi * (6 - math.sqrt(3)) / 23
check("exponent 2π(6−√3)/23", exponent_correct, 1.16593, tol=1e-4)
print(f"          [INFO] exponent = {exponent_correct:.8f}")

scale_factor = math.exp(exponent_correct)
check("exp(2π(6−√3)/23)", scale_factor, 3.2089, tol=1e-3)
print(f"          [INFO] exp(exponent) = {scale_factor:.6f}")

# Wrong exponent from the P127 typo: π(6−√3)/23 (missing factor of 2)
exponent_wrong = math.pi * (6 - math.sqrt(3)) / 23
scale_factor_wrong = math.exp(exponent_wrong)
mu_match_wrong = M_CONF_MEV * scale_factor_wrong
print(f"\n  [INFO] TYPO CHECK: exp(π(6−√3)/23) = {scale_factor_wrong:.6f}")
print(f"          Wrong μ_match = {M_CONF_MEV:.4f} × {scale_factor_wrong:.6f} = {mu_match_wrong:.2f} MeV")
print(f"          (should be ≈ 393–394 MeV, not 705.9 MeV — wrong result confirms typo)")
_mark(_next(),
      f"Wrong exponent gives {mu_match_wrong:.2f} MeV — clearly wrong (differs by "
      f"{abs(mu_match_wrong - 705.9):.1f} MeV from 705.9)",
      abs(mu_match_wrong - 705.9) > 100)

# ────────────────────────────────────────────────────────────────────────────
print("S6  Matching scale μ_match = m_conf · exp(2π(6−√3)/23)")
# ────────────────────────────────────────────────────────────────────────────
MU_MATCH_MEV = M_CONF_MEV * scale_factor
check("μ_match (MeV)", MU_MATCH_MEV, 705.9, tol=1e-3, unit=" MeV")
print(f"          [INFO] μ_match = {MU_MATCH_MEV:.4f} MeV")

# Cross-check the paper's inline arithmetic: 219.99 × 3.209 ≈ 705.9
cross = 219.99 * 3.209
check("219.99 × 3.209 (paper's check)", cross, 705.9, tol=1e-2, unit=" MeV")

# ────────────────────────────────────────────────────────────────────────────
print("S7  Intermediate logarithm (eq. 3 / log-squared form)")
# ────────────────────────────────────────────────────────────────────────────
# ln(μ_match²/m_conf²) = (2√3−1)·12π/(23√3)
rhs_log_sq = (2*math.sqrt(3) - 1) * 12 * math.pi / (23 * math.sqrt(3))
lhs_log_sq = 2 * exponent_correct            # = ln(μ²/μ₀²) = 2 ln(μ/μ₀)
check_exact("ln(μ²/m_conf²): both forms agree", lhs_log_sq, rhs_log_sq, tol=1e-10)
print(f"          [INFO] ln(μ_match²/m_conf²) = {lhs_log_sq:.8f}")

# Algebraic step eq.(4): (2√3−1)/√3 = (6−√3)/3
lhs_alg = (2*math.sqrt(3) - 1) / math.sqrt(3)
rhs_alg = (6 - math.sqrt(3)) / 3
check_exact("(2√3−1)/√3 = (6−√3)/3", lhs_alg, rhs_alg)

# Compact squared-log eq.(5): 4π(6−√3)/23
compact_sq = 4 * math.pi * (6 - math.sqrt(3)) / 23
check_exact("4π(6−√3)/23 = 2×exponent", compact_sq, 2 * exponent_correct)

# ────────────────────────────────────────────────────────────────────────────
print("S8  One-loop RG self-consistency")
# ────────────────────────────────────────────────────────────────────────────
# Starting from α_s(m_conf)=√3, run to μ_match via one-loop MSbar:
#   α_s(μ) = α_s(m_conf) / (1 + b₀/(4π) · α_s(m_conf) · ln(μ²/m_conf²))
log_ratio_sq = 2 * exponent_correct    # ln(μ_match²/m_conf²)
alpha_s_rg = alpha_s_conf / (1 + B0/(4*math.pi) * alpha_s_conf * log_ratio_sq)
check("α_s(μ_match) from one-loop RG", alpha_s_rg, 0.5, tol=1e-6)

# ────────────────────────────────────────────────────────────────────────────
print("S9  Λ_QCD Landau pole (eq. 6)")
# ────────────────────────────────────────────────────────────────────────────
# From m_conf: Λ_QCD = m_conf · exp(−6π/(23√3))
exp_lambda_from_conf = math.exp(-6 * math.pi / (23 * math.sqrt(3)))
LAMBDA_QCD = M_CONF_MEV * exp_lambda_from_conf
check("Λ_QCD from m_conf (MeV)", LAMBDA_QCD, 137.1, tol=1e-2, unit=" MeV")
print(f"          [INFO] Λ_QCD = {LAMBDA_QCD:.4f} MeV")

# From μ_match: Λ_QCD = μ_match · exp(−12π/23)
exp_lambda_from_match = math.exp(-12 * math.pi / 23)
lambda_from_match = MU_MATCH_MEV * exp_lambda_from_match
check("Λ_QCD from μ_match (MeV)", lambda_from_match, 137.1, tol=1e-2, unit=" MeV")

# Both routes give the same Landau pole
check_exact("Λ_QCD both routes agree", LAMBDA_QCD, lambda_from_match, tol=1e-6)

# μ_match / Λ_QCD = exp(12π/23)  (§3 cross-check)
ratio_match_lambda = MU_MATCH_MEV / LAMBDA_QCD
check("μ_match/Λ_QCD = exp(12π/23)", ratio_match_lambda,
      math.exp(12 * math.pi / 23), tol=1e-6)

# ────────────────────────────────────────────────────────────────────────────
print("S10 Factorised exponent forms (eq. 7)")
# ────────────────────────────────────────────────────────────────────────────
# 2π(6−√3)/23  =  2π·ρ_G₂(2ρ_G₂−1) / (3b₀)  =  2π(|Φ⁺|−ρ_G₂) / (3b₀)
form_A = 2 * math.pi * (6 - math.sqrt(3)) / 23
form_B = 2 * math.pi * RHO_G2 * (2*RHO_G2 - 1) / (3 * B0)
form_C = 2 * math.pi * (PHI_POS - RHO_G2) / (3 * B0)
check_exact("form_A = form_B", form_A, form_B)
check_exact("form_A = form_C", form_A, form_C)
print(f"          [INFO] all three exponent forms = {form_A:.10f}")

# ────────────────────────────────────────────────────────────────────────────
print("S11 Golden-ratio near-miss (§4)")
# ────────────────────────────────────────────────────────────────────────────
PHI_GR = (1 + math.sqrt(5)) / 2     # golden ratio
two_phi = 2 * PHI_GR
exp_exact = math.exp(exponent_correct)    # ≈ 3.2089
discrepancy_pct = (two_phi - exp_exact) / exp_exact * 100
check("2φ (golden ratio near-miss)", two_phi, 3.2361, tol=1e-4)
check("exp(exponent) exact", exp_exact, 3.2089, tol=1e-4)
print(f"          [INFO] 2φ = {two_phi:.6f},  exp(exponent) = {exp_exact:.6f}")
print(f"          [INFO] discrepancy = {discrepancy_pct:+.4f}% (paper states +0.85%)")
check("discrepancy 2φ vs exact (%)", discrepancy_pct, 0.85, tol=0.02)
_mark(_next(),
      f"2φ exceeds exact result by {discrepancy_pct:.3f}% — non-structural near-miss confirmed",
      discrepancy_pct > 0)

# ────────────────────────────────────────────────────────────────────────────
print("S12 α_s(M_Z) perturbative predictions (§1, stated results)")
# ────────────────────────────────────────────────────────────────────────────
# These are quoted results from Addendum 127; we check they are in the
# PDG window α_s(M_Z) = 0.1179 ± 0.0009
PDG_AS_MZ   = 0.1179
PDG_AS_ERR  = 0.0009
AS_NLO      = 0.1183   # +0.4σ
AS_N2LO     = 0.1176   # −0.3σ

for label, val, sigma_claim in [("α_s^NLO(M_Z)",   AS_NLO,  +0.4),
                                 ("α_s^N2LO(M_Z)", AS_N2LO, -0.3)]:
    sigma_actual = (val - PDG_AS_MZ) / PDG_AS_ERR
    _mark(_next(),
          f"{label} = {val}  ({sigma_actual:+.1f}σ, paper claims {sigma_claim:+.1f}σ)",
          abs(sigma_actual - sigma_claim) < 0.05)

# Check NLO < N2LO and both within 1σ
assert AS_NLO > PDG_AS_MZ, "NLO should sit above PDG central"
assert AS_N2LO < PDG_AS_MZ, "N2LO should sit below PDG central"
bracket_ok = abs(AS_NLO - PDG_AS_MZ) < PDG_AS_ERR and \
             abs(AS_N2LO - PDG_AS_MZ) < PDG_AS_ERR
_mark(_next(), "NLO / N2LO bracket the PDG central value within 1σ", bracket_ok)

# ════════════════════════════════════════════════════════════════════════════
print(f"\n{'='*60}\nRESULT: {PASS} PASS / {FAIL} FAIL")
sys.exit(0)  # baseline convention: P132 reports but never gates on exit code