"""
verify_P199.py — Numerical verification suite for Addendum 199.

Checks all numerical claims in:
  Lumen/corpus/addenda/199_Addendum_G2_three_loop_N.tex

Run with:  python verify_P199.py
All checks must pass with mpmath at 50-digit precision.

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

import sys

from mpmath import mp, mpf, pi, sqrt, log, zeta, cos, fabs, nstr
mp.dps = 50

# ── Constants ─────────────────────────────────────────────────────────────────

ALPHA_INV = 4*pi**3 + pi**2 + pi
alpha     = 1 / ALPHA_INV
OMEGA_0   = pi**3 / 4

R    = mpf('206.7682830')        # m_mu / m_e  (CODATA target)
tree = mpf('207')                 # integer tree-level

c1 = alpha / (2*pi)              # one-loop  ≡ α/(2π)
c2 = (alpha/pi)**2               # two-loop  ≡ (α/π)²
c3 = (alpha/pi)**3               # three-loop ≡ (α/π)³

val_1loop = tree * (1 - c1)
val_2loop = tree * (1 - c1 + c2)
residual  = val_2loop - R        # negative → formula undershoots

N     = residual / (tree * c3)   # raw formula coefficient (negative)
N_pos = -N                       # positive |N| used in additive correction

# Best G₂/A₂ conjecture (Conjecture 7.1)
N_conj = OMEGA_0 * ALPHA_INV * mpf('11') / 4

# G₂ and A₂ numerical constants
dim_G2 = 14; h_G2  = 4; rank_G2 = 2; Phi_G2 = 12
dim_A2 =  8; h_A2  = 3
C2_7  = mpf('2');    C2_14 = mpf('4')
C2_27 = mpf('16')/3; C2_64 = mpf('10')
T7    = mpf('1')/2;  T14   = mpf('2'); T27 = mpf('9')/2

# QED three-loop coefficient (Laporta–Remiddi)
C3_QED = mpf('197')/144 + pi**2*(2*log(2)-1)/12 - pi**4/360 + 3*zeta(3)/4

# ── Helpers ───────────────────────────────────────────────────────────────────

PASS = FAIL = 0
N_CHECK = 0

def check(desc, cond):
    global PASS, FAIL, N_CHECK
    N_CHECK += 1
    ok = bool(cond)
    PASS += ok
    FAIL += (not ok)
    print(f"  [{'PASS' if ok else 'FAIL'}] {N_CHECK:>2}. {desc}")

def close(desc, a, b, tol):
    """Check |a-b|/max(|b|,1) < tol (same predicate as the original assert)."""
    err = fabs(a - b) / max(fabs(b), mpf('1'))
    check(f"{desc}: rel err = {float(err):.3e} < {tol}", err < tol)

# ══════════════════════════════════════════════════════════════════════════════
# §1 — Fundamental constant checks
# ══════════════════════════════════════════════════════════════════════════════

print("S1  Fundamental constants")

check("C01: 137 < ALPHA_INV < 138", ALPHA_INV > 137 and ALPHA_INV < 138)
close("C02: ALPHA_INV = 137.0363", ALPHA_INV, mpf('137.0363'), 1e-5)
close("C03: ALPHA_INV = 4pi^3 + pi^2 + pi", ALPHA_INV, 4*pi**3 + pi**2 + pi, 1e-40)
check("C04: 0.0072 < alpha < 0.0074", alpha > 0.0072 and alpha < 0.0074)
aopi = alpha/pi
check("C05: 0.00232 < alpha/pi < 0.00234", aopi > 0.00232 and aopi < 0.00234)
close("C06: OMEGA_0 = pi^3/4", OMEGA_0, pi**3/4, 1e-40)
check("C07: 7.75 < OMEGA_0 < 7.76", OMEGA_0 > 7.75 and OMEGA_0 < 7.76)
close("C08: c1 = alpha/(2pi)", c1, alpha/(2*pi), 1e-40)
close("C09: c1 = 0.001161", c1, mpf('0.001161'), 5e-4)
close("C10: c2 = 5.395e-6", c2, mpf('5.395e-6'), 5e-4)
close("C11: c3 = 1.253e-8", c3, mpf('1.253e-8'), 5e-4)

# ══════════════════════════════════════════════════════════════════════════════
# §2 — Two-loop formula and residual (Step A)
# ══════════════════════════════════════════════════════════════════════════════

print("S2  Two-loop formula and residual")

check("C12: val_2loop < R (series undershoots at two loops)", val_2loop < R)
close("C13: val_2loop = 206.7607", val_2loop, mpf('206.7607'), 5e-5)
check("C14: residual < 0 (undershoot)", residual < 0)
close("C15: |residual| = 0.00758", fabs(residual), mpf('0.00758'), 1e-3)
check("C16: N from subtractive formula is negative", N < 0)
close("C17: N_pos near 2921", N_pos, mpf('2921'), 5e-4)
close("C18: N_pos = 2920.838", N_pos, mpf('2920.838'), 1e-5)

# Verify round-trip: formula with N_pos additive gives R
val_check = tree * (1 - c1 + c2 + N_pos*c3)
close("C19: round-trip formula with N_pos gives R", val_check, R, 1e-10)

# ══════════════════════════════════════════════════════════════════════════════
# §3 — P198 one- and two-loop overshoot percentages
# ══════════════════════════════════════════════════════════════════════════════

print("S3  One- and two-loop overshoot percentages")

delta_needed = tree - R                           # = 0.2317170...
one_loop_corr   = tree * c1                       # 207 * α/(2π)
two_loop_corr   = tree * c1 - tree * c2          # net two-loop

overshoot_1 = (one_loop_corr - delta_needed) / delta_needed * 100
overshoot_2 = (two_loop_corr - delta_needed) / delta_needed * 100

close("C20: one-loop overshoot = 3.75%", overshoot_1, mpf('3.75'), 1e-2)
close("C21: two-loop overshoot = 3.27%", overshoot_2, mpf('3.27'), 1e-2)
check("C22: two-loop reduces overshoot", overshoot_2 < overshoot_1)

# ══════════════════════════════════════════════════════════════════════════════
# §4 — G₂ / A₂ structure constants
# ══════════════════════════════════════════════════════════════════════════════

print("S4  G2 / A2 structure constants")

check("C23: dim(G2) - h_A2 = 11", dim_G2 - h_A2 == 11)
check("C24: h_G2 = 4", h_G2 == 4)
check("C25: dim_G2 = 14", dim_G2 == 14)
check("C26: h_A2 = 3", h_A2 == 3)
check("C27: C2_14 = h_G2 (Casimir C2(adj, G2) = dual Coxeter number)", C2_14 == h_G2)
close("C28: Dynkin index T(7) = 1/2", T7, mpf('1')/2, 1e-40)
close("C29: C2(27) = 16/3", C2_27, mpf('16')/3, 1e-40)

# ══════════════════════════════════════════════════════════════════════════════
# §5 — Conjecture 7.1: N_pos ≈ OMEGA_0 × ALPHA_INV × 11/4
# ══════════════════════════════════════════════════════════════════════════════

print("S5  Conjecture 7.1: N_pos = OMEGA_0 * ALPHA_INV * 11/4")

close("C30: N_conj = OMEGA_0 * ALPHA_INV * 11/4", N_conj, OMEGA_0 * ALPHA_INV * 11/4, 1e-40)
check("C31: 2920 < N_conj < 2922", N_conj > 2920 and N_conj < 2922)

rel_err_conj = fabs(N_conj - N_pos) / N_pos
check(f"C32: conjecture rel_err = {float(rel_err_conj):.4e} < 0.02%", rel_err_conj < mpf('0.0002'))

ratio = N_pos / (ALPHA_INV * OMEGA_0)
close("C33: N_pos / (ALPHA_INV * OMEGA_0) = 11/4", ratio, mpf('11')/4, 2e-3)

val_conj_formula = tree * (1 - c1 + c2 + N_conj*c3)
residual_conj = fabs(val_conj_formula - R)
check(f"C34: formula residual with conjecture = {float(residual_conj):.2e} < 2e-6", residual_conj < mpf('2e-6'))
check("C35: conj residual << 2-loop residual", residual_conj < fabs(residual) / 1000)

# ══════════════════════════════════════════════════════════════════════════════
# §6 — OMEGA_0 family (P198 candidates) are far from N_pos
# ══════════════════════════════════════════════════════════════════════════════

print("S6  OMEGA_0 family candidates are far from N_pos")

gap_OMEGA0 = fabs(OMEGA_0 - N_pos) / N_pos
check("C36: OMEGA_0 alone is far from N_pos (>99% off)", gap_OMEGA0 > mpf('0.99'))

pred_k2 = OMEGA_0 * (1 + 2*alpha/pi)
gap_k2 = fabs(pred_k2 - N_pos) / N_pos
check("C37: OMEGA_0*(1+2alpha/pi) far from N_pos (>99% off)", gap_k2 > mpf('0.99'))

# ══════════════════════════════════════════════════════════════════════════════
# §7 — QED comparison
# ══════════════════════════════════════════════════════════════════════════════

print("S7  QED comparison")

check("C38: 2.31 < C3_QED < 2.32 (Laporta-Remiddi)", C3_QED > 2.31 and C3_QED < 2.32)
check("C39: N_pos / C3_QED > 1000 (very different scale)", N_pos / C3_QED > 1000)

ratio_AI = N_pos / ALPHA_INV
close("C40: N_pos / ALPHA_INV = 21.31 (not a simple G2 integer)", ratio_AI, mpf('21.31'), 1e-3)

pi10_32 = pi**10 / 32
rel_pi10 = fabs(pi10_32 - N_pos) / N_pos
check("C41: pi^10/32 within 0.3% of N_pos (second-best candidate)", rel_pi10 < mpf('0.003'))
check("C42: conjecture beats pi^10/32", rel_err_conj < rel_pi10)

# ══════════════════════════════════════════════════════════════════════════════
# §8 — Additional consistency checks
# ══════════════════════════════════════════════════════════════════════════════

print("S8  Additional consistency checks")

prod = OMEGA_0 * ALPHA_INV
check("C43: 1062 < OMEGA_0*ALPHA_INV < 1063", prod > 1062 and prod < 1063)
close("C44: 11/h_G2 = 11/4 (the ratio is purely algebraic)", mpf(dim_G2 - h_A2) / h_G2, mpf('11')/4, 1e-40)
check("C45: val_1loop < val_2loop (two-loop term raises the value)", val_1loop < val_2loop)
check("C46: val_2loop < R (two-loop still undershoots)", val_2loop < R)
check("C47: val_conj_formula > R (conjecture overshoots marginally)", val_conj_formula > R)

three_loop_correction = N_pos * c3
check("C48: 1e-5 < three-loop correction < 1e-3",
      three_loop_correction > mpf('1e-5') and three_loop_correction < mpf('1e-3'))
check("C49: c1 > c2 > c3 (loop expansion is ordered)", c1 > c2 > c3)
check("C50: |Phi(G2)| = 12 and rank(G2) = 2", Phi_G2 == 12 and rank_G2 == 2)

# ── Key values (informational) ────────────────────────────────────────────────

print(f"  N_pos  = {nstr(N_pos, 10)}")
print(f"  N_conj = {nstr(N_conj, 10)}  [OMEGA_0 * ALPHA_INV * 11/4]")
print(f"  rel err = {float(fabs(N_conj-N_pos)/N_pos)*100:.5f}%")
print(f"  formula residual (conj) = {float(fabs(val_conj_formula-R)):.3e}")

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