#!/usr/bin/env python3
"""
verify_P124.py
Computational verification of all numerical claims in:

  Addendum P124 — W-Boson NLO Programme:
  mt Scheme Resolution and Electroweak Corrections Beyond Delta-rho.

Every value is re-derived from first principles using TOE constants.
Claimed values are compared with independently computed values; each
check is annotated with absolute error, relative error, and PASS/FAIL.

Author  : Leon Fernando Vlegels
License : MIT
"""

import math
import sys

# ===========================================================================
# TOE mathematical constants
# ===========================================================================

ALPHA_INV     = 4 * math.pi**3 + math.pi**2 + math.pi   # ≈ 137.036
ALPHA         = 1 / ALPHA_INV                             # fine-structure constant
BREATH_PERIOD = math.pi * ALPHA_INV                       # ≈ 432 system units
PHI           = (1 + math.sqrt(5)) / 2                    # golden ratio ≈ 1.618034

# Weinberg angle from J3(O) G2×2I structure (P112)
SIN2_THETA_W  = 3 / (8 * PHI)   # ≈ 0.23173 from formula

# ===========================================================================
# PDG 2024 physical inputs  (P124 §2 ground-truth table)
# ===========================================================================

MW_PDG   = 80.377          # GeV   (PDG 2024, ±0.012 GeV)
MZ       = 91.1876         # GeV   (PDG 2024)
GF       = 1.16637e-5      # GeV⁻² (PDG 2024; paper uses 1.16637×10⁻⁵ throughout)
MT_PDG   = 172.69          # GeV   (on-shell pole mass, PDG 2024, ±0.30 GeV)
MT_TOE   = 176.101         # GeV   (TOE spectral / GUT-scale running mass, P100)
ALPHA_S  = 0.108           # alpha_s(mt)  (PDG / P106)

# Note: task prompt lists GF = 1.1663788e-5 (higher precision PDG value).
# P124 rounds this to 1.16637e-5; both agree to all digits shown in the paper.

# ===========================================================================
# Verification infrastructure
# ===========================================================================

_results: list = []


def check(label: str,
          computed: float,
          claimed: float,
          tol_rel: float = 5e-3) -> bool:
    """
    Assert |computed - claimed| / |claimed| <= tol_rel.
    Prints a PASS/FAIL line and accumulates results.
    """
    if claimed == 0:
        err_abs = computed - claimed
        err_rel = abs(err_abs)
    else:
        err_abs = computed - claimed
        err_rel = err_abs / abs(claimed)

    ok = abs(err_rel) <= tol_rel
    _results.append((label, ok, computed, claimed, err_rel))
    print(f"  [{'PASS' if ok else 'FAIL'}] {len(_results):>2}. {label}")
    print(f"       computed = {computed:.9g}  |  claimed = {claimed:.9g}"
          f"  |  err = {err_abs:+.3e}  ({err_rel*100:+.4f}%)")
    return ok


def note(text: str) -> None:
    print(f"  [i] {text}")


# ===========================================================================
# Section 0 — Constant cross-checks
# ===========================================================================

print()
print("=" * 72)
print("Section 0 — TOE constants and setup")
print("=" * 72)

note(f"ALPHA_INV     = {ALPHA_INV:.8f}   (expect ≈ 137.036)")
note(f"ALPHA         = {ALPHA:.8e}")
note(f"BREATH_PERIOD = {BREATH_PERIOD:.6f}   (expect ≈ 432)")
note(f"PHI           = {PHI:.12f}")
note(f"SIN2_THETA_W  = 3/(8φ) = {SIN2_THETA_W:.8f}")

# ------ sin²θW table-value discrepancy ------
SIN2_TW_TABLE = 0.23099   # value quoted in P124 §2 ground-truth table
discrepancy   = SIN2_THETA_W - SIN2_TW_TABLE
mw_tree_from_table = MZ * math.sqrt(1 - SIN2_TW_TABLE)

print()
print("  --- sin²θW discrepancy report ---")
print(f"  Formula 3/(8φ):            {SIN2_THETA_W:.8f}")
print(f"  P124 table entry:          {SIN2_TW_TABLE:.8f}")
print(f"  Difference (formula−table): {discrepancy:+.8f}  "
      f"({discrepancy/SIN2_TW_TABLE*100:+.4f}%)")
print(f"  MW_tree using table value:  {mw_tree_from_table:.4f} GeV  ← does NOT match 79.925")
MW_tree = MZ * math.sqrt(1 - SIN2_THETA_W)
print(f"  MW_tree using formula:      {MW_tree:.4f} GeV  ← matches 79.925")
print("  CONCLUSION: the table entry '0.23099' is a numerical error in P124.")
print("  All subsequent calculations use the formula value 3/(8φ).")


# ===========================================================================
# Section 1 — Tree-level MW and P122 baseline
# ===========================================================================

print()
print("=" * 72)
print("Section 1 — Tree-level MW and P122 Δρ baseline")
print("=" * 72)

# (1a) MW_tree = MZ · √(1 − sin²θW)   eq.(1)
MW_tree = MZ * math.sqrt(1 - SIN2_THETA_W)
check("MW_tree = MZ·√(1 − 3/(8φ))", MW_tree, 79.925, tol_rel=2e-3)

# (1b) Δρ_TOE with mt = mt_TOE = 176.101 GeV   eqs.(2–3)
drho_toe_num = 3 * GF * MT_TOE**2
drho_toe_den = 8 * math.pi**2 * math.sqrt(2)
drho_toe     = drho_toe_num / drho_toe_den

# P124 eqs.(2–3) do not quote the TOE numerator explicitly; only the final ratio.
# (1.04353 is the PDG-mass numerator from Theorem 2.1 — checked in Section 2.)
note(f"3·GF·mt_TOE² (numerator, not quoted in paper) = {drho_toe_num:.7f}")
check("denominator 8π²√2",        drho_toe_den, 111.653, tol_rel=1e-4)
check("Δρ_TOE",                   drho_toe,     0.009718, tol_rel=5e-3)

# (1c) MW_P122 = MW_tree · √(1 + Δρ_TOE)   eq.(4)
MW_P122 = MW_tree * math.sqrt(1 + drho_toe)
check("MW_P122 = MW_tree·√(1 + Δρ_TOE)", MW_P122, 80.313, tol_rel=2e-3)

# (1d) P122 residual
res_P122_abs = MW_P122 - MW_PDG
res_P122_pct = res_P122_abs / MW_PDG * 100
check("P122 residual  [GeV]",  res_P122_abs, -0.064, tol_rel=5e-2)
check("P122 residual  [%]",    res_P122_pct, -0.080, tol_rel=5e-2)

# (1e) signed gap entering P124   eq.(5)
delta_0 = MW_P122 - MW_PDG
check("δ₀ = MW_P122 − MW_PDG  [GeV]", delta_0, -0.064, tol_rel=5e-2)


# ===========================================================================
# Section 2 — Step 1: scheme resolution with PDG pole mass (Theorems 2.1–2.2)
# ===========================================================================

print()
print("=" * 72)
print("Section 2 — Step 1: scheme resolution with PDG pole mass")
print("=" * 72)

# (2a) mt_PDG²   eq. block of Theorem 2.1
mt2_pdg = MT_PDG**2
check("mt_PDG² = 172.69²  [GeV²]", mt2_pdg, 29_821.84, tol_rel=1e-5)

# (2b) Numerator and denominator of Δρ_PDG
drho_pdg_num = 3 * GF * mt2_pdg
drho_pdg_den = 8 * math.pi**2 * math.sqrt(2)
check("3·GF·mt_PDG²  (numerator)",   drho_pdg_num, 1.04353, tol_rel=5e-3)
check("8π²√2  (denominator)",        drho_pdg_den, 111.653, tol_rel=1e-4)

# sub-checks on individual factors the paper quotes
check("8 × π² = 8 × 9.86960",    8 * math.pi**2, 8 * 9.86960, tol_rel=1e-5)
check("√2 = 1.41421",             math.sqrt(2),   1.41421,     tol_rel=1e-5)

# (2c) Δρ_PDG   eq.(6)
drho_pdg = drho_pdg_num / drho_pdg_den
check("Δρ_PDG = 3·GF·mt_PDG²/(8π²√2)", drho_pdg, 0.009347, tol_rel=5e-3)

# (2d) MW_LO = MW_tree · √(1 + Δρ_PDG)   Theorem 2.2 / eq.(7)
sqrt_lo = math.sqrt(1 + drho_pdg)
check("√(1 + Δρ_PDG)",                          sqrt_lo, 1.004663, tol_rel=1e-4)
MW_LO = MW_tree * sqrt_lo
check("MW_LO = MW_tree·√(1 + Δρ_PDG)  [GeV]",  MW_LO,   80.298,   tol_rel=2e-3)

# (2e) LO residual
res_LO_abs = MW_LO - MW_PDG
res_LO_pct = res_LO_abs / MW_PDG * 100
check("LO residual  [GeV]",  res_LO_abs, -0.079, tol_rel=5e-2)
check("LO residual  [%]",    res_LO_pct, -0.099, tol_rel=5e-2)

# (2f) Scheme-correction opening the gap (P122→P124 LO)
scheme_delta = MW_LO - MW_P122
note(f"Scheme correction MW_LO − MW_P122 = {scheme_delta*1000:.1f} MeV  "
     f"(paper says ~−15 MeV)")

# (2g) Informational: mt running estimate (not a precision claim)
print()
note("-- mt running estimate (informational, ±2% accuracy) --")
mt_msbar_est = MT_TOE * 0.970
note(f"mt_MS̄(mt) ≈ mt_TOE × 0.970 = {mt_msbar_est:.1f} GeV  (claimed ≈ 170.8)")
pole_factor_corr = 4 * ALPHA_S / (3 * math.pi)
check("4·αs/(3π)  (one-loop pole-mass factor)",
      pole_factor_corr, 0.04584, tol_rel=5e-3)
check("1 + 4·αs/(3π)  (paper quotes 1.0459)",
      1 + pole_factor_corr, 1.0459, tol_rel=5e-3)
mt_pole_est = mt_msbar_est * (1 + pole_factor_corr)
check("mt_pole estimate = mt_MS̄·(1 + 4αs/(3π))  [GeV]",
      mt_pole_est, 178.6, tol_rel=5e-3)


# ===========================================================================
# Section 3 — QCD correction to Δρ at O(αs)   (§3.1)
# ===========================================================================

print()
print("=" * 72)
print("Section 3 — QCD correction to Δρ at O(αs)")
print("=" * 72)

# (3a) QCD factor 8·αs/(3π)
qcd_factor = 8 * ALPHA_S / (3 * math.pi)
check("8·αs(mt)/(3π)",           qcd_factor, 0.09166, tol_rel=2e-3)
note(f"Numerator  8×0.108  = {8*ALPHA_S:.4f}  (paper: 0.864)")
note(f"Denominator  3π     = {3*math.pi:.5f}  (paper: 9.4248)")
check("numerator  8×αs",     8 * ALPHA_S,    0.864,  tol_rel=1e-9)
check("denominator 3π",      3 * math.pi,    9.4248, tol_rel=1e-4)

# (3b) Δρ_QCD = Δρ_PDG · (1 − 8αs/(3π))   eq.(9)
drho_qcd   = drho_pdg * (1 - qcd_factor)
one_m_qcd  = 1 - qcd_factor
check("(1 − 8αs/(3π))",                        one_m_qcd, 0.90834, tol_rel=2e-3)
check("Δρ_QCD = Δρ_PDG·(1 − 8αs/(3π))",       drho_qcd,  0.008490, tol_rel=2e-3)

# (3c) δΔρ_QCD   eq.(10)
d_drho_qcd = drho_qcd - drho_pdg
check("δΔρ_QCD = Δρ_QCD − Δρ_PDG",            d_drho_qcd, -0.000857, tol_rel=5e-3)

# (3d) δMW_QCD ≈ (MW_tree/2) · δΔρ_QCD   eq.(11)  [linearised]
dMW_qcd = (MW_tree / 2) * d_drho_qcd
check("δMW_QCD  [GeV]   (linearised)",         dMW_qcd,         -0.034, tol_rel=5e-2)
check("δMW_QCD  [MeV]   (linearised)",         dMW_qcd * 1000,  -34,    tol_rel=5e-2)

# (3e) MW_QCD = MW_tree · √(1 + Δρ_QCD)   (full formula, eq. after (11))
sqrt_qcd = math.sqrt(1 + drho_qcd)
check("√(1 + Δρ_QCD)",                         sqrt_qcd, 1.004237, tol_rel=1e-4)
MW_QCD = MW_tree * sqrt_qcd
check("MW_QCD = MW_tree·√(1 + Δρ_QCD)  [GeV]", MW_QCD,   80.264,   tol_rel=2e-3)

# (3f) Residual after QCD correction
res_QCD_abs = MW_QCD - MW_PDG
res_QCD_pct = res_QCD_abs / MW_PDG * 100
check("Residual after QCD  [GeV]",  res_QCD_abs, -0.113, tol_rel=5e-2)
check("Residual after QCD  [%]",    res_QCD_pct, -0.141, tol_rel=5e-2)


# ===========================================================================
# Section 4 — Subleading oblique at O(GF² mt⁴)   (§3.2)
# ===========================================================================

print()
print("=" * 72)
print("Section 4 — Subleading oblique correction O(GF² mt⁴)")
print("=" * 72)

# (4a) δΔρ⁽²⁾ ≈ (Δρ_PDG)²   eq.(12)
d_drho_subl = drho_pdg**2
check("δΔρ⁽²⁾ = (Δρ_PDG)²",          d_drho_subl, 8.74e-5, tol_rel=5e-3)

# (4b) δMW⁽²⁾ ≈ (MW_tree/2) · δΔρ⁽²⁾   eq.(13)
dMW_subl = (MW_tree / 2) * d_drho_subl
check("δMW⁽²⁾  [GeV]",                dMW_subl,        0.0035, tol_rel=5e-2)
check("δMW⁽²⁾  [MeV]",                dMW_subl * 1000,  4,     tol_rel=2e-1)

# (4c) MW after QCD + subleading (gap-closure table row 4)
MW_qcd_subl = MW_LO + dMW_qcd + dMW_subl
check("MW (QCD + subleading oblique)  [GeV]",  MW_qcd_subl, 80.268, tol_rel=2e-3)
res_qs_abs = MW_qcd_subl - MW_PDG
res_qs_pct = res_qs_abs / MW_PDG * 100
check("Residual QCD + subleading  [GeV]",  res_qs_abs, -0.109, tol_rel=5e-2)
check("Residual QCD + subleading  [%]",    res_qs_pct, -0.136, tol_rel=5e-2)


# ===========================================================================
# Section 5 — EW vertex+box corrections at O(α/π)   (§3.3)
# ===========================================================================

print()
print("=" * 72)
print("Section 5 — EW vertex+box corrections O(α/π)  [SM literature input]")
print("=" * 72)

# This correction is taken directly from SM calculations (Sirlin 1980,
# Denner 1991) and cannot be independently re-derived from TOE constants
# alone.  The verification below checks only the arithmetic use of the
# stated value.
dMW_ew = -0.038  # GeV  (−38 MeV, from SM literature as cited in P124)
note("δMW_EW = −38 MeV  (SM lit. input: Sirlin 1980 / Denner 1991)")
note("This value is accepted as given; only its downstream arithmetic is verified.")
check("δMW_EW  [MeV]  (identity check)",  dMW_ew * 1000, -38, tol_rel=1e-9)


# ===========================================================================
# Section 6 — NLO sum and MW_NLO   (§4 / Theorem 4.1)
# ===========================================================================

print()
print("=" * 72)
print("Section 6 — NLO sum and final MW_NLO  (Theorem 4.1)")
print("=" * 72)

# (6a) integer MeV arithmetic   eq.(14)
check("−34 + 4 + (−38) = −68  [MeV]",   -34 + 4 + (-38), -68, tol_rel=1e-9)

# (6b) sum of computed shifts
dMW_nlo = dMW_qcd + dMW_subl + dMW_ew
check("δMW_NLO = δMW_QCD + δMW⁽²⁾ + δMW_EW  [GeV]", dMW_nlo,       -0.068, tol_rel=5e-2)
check("δMW_NLO  [MeV]",                               dMW_nlo * 1000, -68,   tol_rel=5e-2)

# (6c) individual MeV components
check("δMW_QCD  [MeV]",  dMW_qcd * 1000,  -34, tol_rel=5e-2)
check("δMW⁽²⁾  [MeV]",   dMW_subl * 1000,  +4, tol_rel=2e-1)
check("δMW_EW   [MeV]",   dMW_ew * 1000,  -38, tol_rel=1e-9)

# (6d) MW_NLO = MW_LO + δMW_NLO   eq.(15)
MW_NLO = MW_LO + dMW_nlo
check("MW_NLO = MW_LO + δMW_NLO  [GeV]",  MW_NLO, 80.230, tol_rel=2e-3)

# (6e) NLO residual   eq.(16)
res_NLO_abs = MW_NLO - MW_PDG
res_NLO_pct = res_NLO_abs / MW_PDG * 100
check("NLO residual  [GeV]",  res_NLO_abs, -0.147, tol_rel=5e-2)
check("NLO residual  [%]",    res_NLO_pct, -0.183, tol_rel=5e-2)

# (6f) Explicit arithmetic from eq.(15): 80.298 + (−0.034) + (+0.004) + (−0.038)
arithmetic_check = 80.298 + (-0.034) + 0.004 + (-0.038)
check("80.298 − 0.034 + 0.004 − 0.038 = 80.230",
      arithmetic_check, 80.230, tol_rel=1e-5)


# ===========================================================================
# Section 7 — Gap-closure table (§4.2)
# ===========================================================================

print()
print("=" * 72)
print("Section 7 — Gap-closure table: all five rows")
print("=" * 72)

rows = [
    ("Row 1 — P122 (TOE mt, Δρ only)",      MW_P122,      80.313, -0.064, -0.080),
    ("Row 2 — P124 LO (PDG mt, Δρ only)",   MW_LO,        80.298, -0.079, -0.099),
    ("Row 3 — After QCD correction to Δρ",  MW_QCD,       80.264, -0.113, -0.141),
    ("Row 4 — After QCD + subleading",       MW_qcd_subl,  80.268, -0.109, -0.136),
    ("Row 5 — After all NLO  (Theorem 4.1)", MW_NLO,       80.230, -0.147, -0.183),
]

for label, MW_comp, mw_claim, res_claim_gev, res_claim_pct in rows:
    print(f"\n  {label}:")
    check(f"  MW  [GeV]",  MW_comp,             mw_claim,       tol_rel=2e-3)
    check(f"  residual [GeV]", MW_comp - MW_PDG, res_claim_gev,  tol_rel=5e-2)
    check(f"  residual [%]",
          (MW_comp - MW_PDG) / MW_PDG * 100, res_claim_pct,    tol_rel=5e-2)


# ===========================================================================
# Section 8 — sin²θW numerical value (table cross-check)
# ===========================================================================

print()
print("=" * 72)
print("Section 8 — sin²θW formula vs. table-value cross-check")
print("=" * 72)

note(f"sin²θW = 3/(8φ)       = {SIN2_THETA_W:.8f}   (formula)")
note(f"sin²θW  P124 table    = {SIN2_TW_TABLE:.8f}   (stated in §2 table)")
note(f"Difference             = {SIN2_THETA_W - SIN2_TW_TABLE:+.8f}")
note(f"MW_tree from table val = {mw_tree_from_table:.4f} GeV  (≠ 79.925 GeV)")
note(f"MW_tree from formula   = {MW_tree:.4f} GeV  (≈ 79.925 GeV ✓)")
note("The table entry 0.23099 is numerically inconsistent with:")
note("  (a) the formula 3/(8φ), and")
note("  (b) the subsequent computation MW_tree = 79.925 GeV.")
note("All arithmetic in this verifier uses the formula value.")

# Verify the formula sin²θW = 3/(8φ) = 3(√5−1)/16  (algebraic identity)
sin2_algebraic = 3 * (math.sqrt(5) - 1) / 16
check("3/(8φ) = 3(√5−1)/16  (algebraic identity)",
      SIN2_THETA_W, sin2_algebraic, tol_rel=1e-12)


# ===========================================================================
# Final summary
# ===========================================================================

print()
print("=" * 72)
print("FINAL SUMMARY")
print("=" * 72)

passes = sum(1 for _, ok, *_ in _results if ok)
fails  = sum(1 for _, ok, *_ in _results if not ok)
total  = len(_results)

print()
print(f"  Arithmetic checks: {passes} PASS / {fails} FAIL  (of {total} total)")

if fails:
    print()
    print("  FAILED checks:")
    for label, ok, comp, claim, rel_err in _results:
        if not ok:
            print(f"    ✗  {label}")
            print(f"       computed={comp:.9g}, claimed={claim:.9g}, "
                  f"rel_err={rel_err*100:+.4f}%")

print()
print("  Key derived values:")
print(f"    PHI                        = {PHI:.10f}")
print(f"    sin²θW = 3/(8φ)            = {SIN2_THETA_W:.8f}  "
      f"(table claims 0.23099 — see §8)")
print(f"    MW_tree                    = {MW_tree:.4f} GeV  (claimed 79.925)")
print(f"    Δρ_TOE  [P122, mt=176.101] = {drho_toe:.7f}     (claimed 0.009718)")
print(f"    MW_P122                    = {MW_P122:.4f} GeV  (claimed 80.313)")
print(f"    Δρ_PDG  [PDG pole mass]    = {drho_pdg:.7f}     (claimed 0.009347)")
print(f"    MW_LO                      = {MW_LO:.4f} GeV  (claimed 80.298)")
print(f"    δMW_QCD                    = {dMW_qcd*1000:.2f} MeV  (claimed −34 MeV)")
print(f"    δMW_subl                   = {dMW_subl*1000:.2f} MeV  (claimed +4 MeV)")
print(f"    δMW_EW                     = {dMW_ew*1000:.1f} MeV  (SM lit. input)")
print(f"    δMW_NLO                    = {dMW_nlo*1000:.2f} MeV  (claimed −68 MeV)")
print(f"    MW_NLO                     = {MW_NLO:.4f} GeV  (claimed 80.230)")
print(f"    NLO residual               = {res_NLO_abs:.4f} GeV  (claimed −0.147)")
print(f"    NLO residual               = {res_NLO_pct:.4f}%  (claimed −0.183%)")
print()
print("  Caveats:")
print("    • δMW_EW = −38 MeV is a SM literature input (Sirlin 1980 / Denner 1991).")
print("      Its sign and magnitude are verified arithmetically but not re-derived")
print("      from TOE constants, as P124 explicitly defers this to B-1/B-2/B-3.")
print("    • sin²θW table entry 0.23099 is a numerical error in P124 §2;")
print("      the formula 3/(8φ) = 0.23173 is consistent with all MW computations.")
print()

if fails == 0:
    print("  OVERALL: PASS — all arithmetic claims verified within tolerance.")
else:
    print(f"  OVERALL: FAIL — {fails} check(s) failed.")

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