"""
verify_P198.py — Addendum 198: G2 loop integral and mass formula correction.

Checks all key computations from Steps A-F: Delta, alpha/2pi, ratio r,
candidate corrections, two-loop formula, root-angle candidates, and
reverse-engineering identities.  Requires mpmath.

Usage:  python verify_P198.py
All 30+ assert statements must pass for the addendum to be considered
numerically consistent.
"""

import math
try:
    import mpmath as mp
    mp.mp.dps = 50
    HAS_MPMATH = True
except ImportError:
    HAS_MPMATH = False
    import math as mp   # fallback — reduced precision


# ===========================================================
# Constants (TOE exact values)
# ===========================================================
if HAS_MPMATH:
    ALPHA_INV = 4*mp.pi**3 + mp.pi**2 + mp.pi
    alpha     = mp.mpf(1) / ALPHA_INV
    OMEGA_0   = mp.pi**3 / 4
    BREATH_PERIOD = mp.pi * ALPHA_INV
    sqrt3     = mp.sqrt(3)
    pi        = mp.pi
else:
    ALPHA_INV = 4*math.pi**3 + math.pi**2 + math.pi
    alpha     = 1.0 / ALPHA_INV
    OMEGA_0   = math.pi**3 / 4
    BREATH_PERIOD = math.pi * ALPHA_INV
    sqrt3     = math.sqrt(3)
    pi        = math.pi

# G2 root-system data
hv_G2   = 4      # dual Coxeter number
dim_G2  = 14     # dimension
roots_G2 = 12    # number of roots
hv_A2   = 3      # dual Coxeter of A2 ⊂ G2
dim_A2  = 8      # dimension of A2
roots_A2 = 6     # number of roots of A2
C2_fund_G2 = 2   # Casimir of 7-dim fund. rep
C2_adj_G2  = 4   # Casimir of adjoint = h^v

# PDG mass ratio and tree-level formula
m_ratio = float(206.7682830)
tree    = float(207)               # dim(G2)*(dim(G2)+1) - dim(A2) = 14*15 - 3

tol_rel = 1e-7   # relative tolerance for float comparisons
tol_abs = 1e-8   # absolute tolerance

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

def rel_err(a, b):
    return abs(float(a) - float(b)) / max(abs(float(b)), 1e-300)

def abs_err(a, b):
    return abs(float(a) - float(b))


# ===========================================================
# §0 — Structural checks
# ===========================================================
# Check 1: tree-level formula
check(1, "tree-level formula: 207 = 14*15 - 3", tree == 14*15 - 3)

# Check 2: ALPHA_INV magnitude
check(2, "ALPHA_INV magnitude in (137.035, 137.037)", 137.035 < float(ALPHA_INV) < 137.037)

# Check 3: alpha = 1/ALPHA_INV
check(3, "alpha = 1/ALPHA_INV roundtrip", rel_err(alpha, 1.0/float(ALPHA_INV)) < tol_rel)

# Check 4: OMEGA_0 = pi^3/4
check(4, "OMEGA_0 = pi^3/4", rel_err(OMEGA_0, math.pi**3/4) < tol_rel)

# Check 5: BREATH_PERIOD = pi * ALPHA_INV
check(5, "BREATH_PERIOD = pi * ALPHA_INV", rel_err(BREATH_PERIOD, math.pi * float(ALPHA_INV)) < tol_rel)

# Check 6: G2 Casimir consistency: C2(adj) = h^v
check(6, "G2 Casimir consistency: C2(adj) = h^v", C2_adj_G2 == hv_G2)

# Check 7: dim_G2 = 14, roots_G2 = 12
check(7, "dim_G2 = 14, roots_G2 = 12", dim_G2 == 14 and roots_G2 == 12)

# Check 8: A2 embedding: roots_A2 = 6, dim_A2 = 8
check(8, "A2 embedding: roots_A2 = 6, dim_A2 = 8", roots_A2 == 6 and dim_A2 == 8)

# Check 9: long/short ratio squared = 3
check(9, "long/short ratio squared = 3", abs(float(sqrt3)**2 - 3) < 1e-12)

# Check 10: m_ratio is in expected PDG range
check(10, "m_ratio in expected PDG range", 206.766 < m_ratio < 206.770)


# ===========================================================
# §A — Delta
# ===========================================================
Delta = 1.0 - m_ratio / tree

# Check 11: Delta positive and small
check(11, "Delta positive and small (1.0e-3, 1.5e-3)", 1.0e-3 < Delta < 1.5e-3)

# Check 12: 207 * Delta = 207 - m_ratio
check(12, "207 * Delta = 207 - m_ratio", abs_err(tree * Delta, tree - m_ratio) < 1e-10)

# Check 13: Delta to 10 significant figures
check(13, "Delta = 1.119405797e-3 to 10 significant figures", rel_err(Delta, 1.119405797e-3) < 1e-6)

# Check 14: 207 * Delta numerically
check(14, "207 * Delta = 0.2317170 numerically", rel_err(tree * Delta, 0.2317170) < 1e-5)

# alpha/(2pi)
a_2pi = float(alpha) / (2 * math.pi)

# Check 15: alpha/(2pi) in expected range
check(15, "alpha/(2pi) in expected range", 1.15e-3 < a_2pi < 1.18e-3)

# Check 16: alpha/(2pi) precise value
check(16, "alpha/(2pi) = 1.16140715e-3 precise value", rel_err(a_2pi, 1.16140715e-3) < 1e-6)

# Overshoot ratio r = alpha/(2pi) / Delta
r = a_2pi / Delta

# Check 17: r > 1 (one-loop overshoots)
check(17, "r > 1 (one-loop overshoots)", r > 1.0)

# Check 18: r precise value
check(18, "r = 1.037521 precise value", rel_err(r, 1.037521) < 1e-4)

# Check 19: r - 1 ~ 3.75%
check(19, "r - 1 ~ 3.75%", 0.036 < (r - 1) < 0.040)


# ===========================================================
# §B — Candidate corrections
# ===========================================================

def gap(c):
    """Return |207*(1-c) - m_ratio|."""
    return abs(tree * (1 - c) - m_ratio)

def predicted(c):
    return tree * (1 - c)

# --- Candidate 1: alpha/(2pi) ---
c1 = a_2pi
# Check 20: c1 gap < 0.01
check(20, "c1 = alpha/(2pi): gap < 0.01", gap(c1) < 0.01)

# Check 21: c1 overshoots Delta
check(21, "c1 overshoots Delta", c1 > Delta)

# --- Candidate 2: alpha/(2pi)*(1-alpha/pi) ---
c2 = a_2pi * (1 - float(alpha)/math.pi)
# Check 22: c2 < c1 (one-loop sub-correction)
check(22, "c2 < c1 (one-loop sub-correction)", c2 < c1)

# Check 23: c2 gap < c1 gap
check(23, "c2 gap < c1 gap (two-step correction improves)", gap(c2) < gap(c1))

# --- Candidate 3 (best): alpha/(2pi)*(1-2alpha/pi) = two-loop QED ---
c3 = a_2pi * (1 - 2*float(alpha)/math.pi)
c3_alt = a_2pi - (float(alpha)/math.pi)**2  # algebraically identical

# Check 24: algebraic identity c3 == c3_alt
check(24, "algebraic identity c3 == c3_alt", abs_err(c3, c3_alt) < 1e-15)

# Check 25: c3 is best Step B candidate (smallest gap)
check(25, "c3 is best Step B candidate (beats c1 and c2)", gap(c3) < gap(c1) and gap(c3) < gap(c2))

# Check 26: c3 gap precise value
check(26, "c3 gap = 7.577e-3 precise value", rel_err(gap(c3), 7.577e-3) < 0.01)

# Check 27: c3 predicted value
check(27, "c3 predicted value = 206.76071", rel_err(predicted(c3), 206.76071) < 1e-4)

# --- G2-specific candidates ---
c_2alpha_7pi = 2*float(alpha)/(7*math.pi)
c_alpha_4pi  = float(alpha)/(4*math.pi)
c_3alpha_14pi = 3*float(alpha)/(14*math.pi)

# Check 28: these G2-weighted candidates all worse than c3
check(28, "G2-weighted candidates (2a/7pi, a/4pi, 3a/14pi) all worse than c3",
      gap(c_2alpha_7pi) > gap(c3) and gap(c_alpha_4pi) > gap(c3)
      and gap(c_3alpha_14pi) > gap(c3))

# Root-angle candidate: alpha*cos(pi/6)/(2pi)
c_cos30 = float(alpha)*math.cos(math.pi/6)/(2*math.pi)
# Check 29: cos(pi/6) candidate meaningfully worse than c3 (>2x)
check(29, "cos(pi/6) candidate meaningfully worse than c3 (>2x)", gap(c_cos30) > 2*gap(c3))

# 1/BREATH_PERIOD candidate
c_breath = 1.0 / float(BREATH_PERIOD)
# Check 30: BREATH_PERIOD-based correction is not competitive
check(30, "BREATH_PERIOD-based correction is not competitive", gap(c_breath) > 0.1)


# ===========================================================
# §C — Two-loop structure
# ===========================================================

# Check 31: two-loop reduces overshoot vs one-loop
check(31, "two-loop reduces overshoot vs one-loop", gap(c3) < gap(c1))

# Check 32: the residual epsilon = Delta - c3 is negative (overshoots)
epsilon = Delta - c3
check(32, "residual epsilon = Delta - c3 is negative (overshoots)", epsilon < 0)

# Check 33: |epsilon| ~ 3.66e-5
check(33, "|epsilon| ~ 3.66e-5", rel_err(abs(epsilon), 3.66e-5) < 0.05)

# Check 34: epsilon is O(alpha^2) — between (a/pi)^2 and (a/pi)
check(34, "epsilon is O(alpha^2) -- between (a/pi)^2 and (a/pi)", (float(alpha)/math.pi)**2 < abs(epsilon) < float(alpha)/math.pi)

# G2 variant: alpha/(2pi)*(1 - 4*alpha/(3pi))
c_G2v2 = a_2pi * (1 - 4*float(alpha)/(3*math.pi))
# Check 35: G2 variant is worse than c3
check(35, "G2 two-loop variant worse than standard c3", gap(c_G2v2) > gap(c3))


# ===========================================================
# §D — Root angle geometry
# ===========================================================

# Check 36: sqrt(3) = tan(pi/3) and G2 root-length ratio
check(36, "sqrt(3) = tan(pi/3) and G2 root-length ratio", abs_err(sqrt3, math.tan(math.pi/3)) < 1e-12)

# Check 37: min angle between G2 roots is pi/6 = 30 deg
check(37, "min angle between G2 roots is pi/6: cos(pi/6) = sqrt(3)/2", abs_err(math.cos(math.pi/6), math.sqrt(3)/2) < 1e-12)

# Check 38: all root-angle corrections worse than c3
c_cos_pi6 = float(alpha)*math.cos(math.pi/6)/(2*math.pi)
c_2_sqrt3  = 2*float(alpha)/(3*math.pi*math.sqrt(3))
check(38, "all root-angle corrections worse than c3",
      gap(c_cos_pi6) > gap(c3) and gap(c_2_sqrt3) > gap(c3))


# ===========================================================
# §E/F — Reverse-engineering
# ===========================================================
D = tree - m_ratio  # absolute deficit

# Check 39: D = 207*Delta
check(39, "D = 207*Delta", rel_err(D, tree * Delta) < 1e-8)

# Check 40: D to 5 sig figs
check(40, "D = 0.231717 to 5 sig figs", rel_err(D, 0.231717) < 1e-4)

# Check 41: D/alpha known dimensionless ratio
Da = D / float(alpha)
check(41, "D/alpha = 31.754 dimensionless ratio", rel_err(Da, 31.754) < 1e-3)

# x_factor
x_factor = 1.0 - Delta * 2 * math.pi / float(alpha)
# Check 42: x_factor > 0 (overshoot)
check(42, "x_factor > 0 (overshoot)", x_factor > 0)

# Check 43: x_factor ~ 0.03616
check(43, "x_factor ~ 0.036164", rel_err(x_factor, 0.036164) < 1e-3)

# N = (a/2pi - Delta)/(a/pi)^2
N = (a_2pi - Delta) / (float(alpha)/math.pi)**2
# Check 44: N computed correctly
check(44, "N = (a/2pi - Delta)/(a/pi)^2 = 7.784", rel_err(N, 7.784) < 1e-2)

# Check 45: N is not a simple G2 Casimir (not 4, 8, 12, 14)
G2_integers = [1, 2, 3, 4, 6, 7, 8, 12, 14]
check(45, "N is not a simple G2 Casimir (not 4, 8, 12, 14)", all(abs(N - k) > 0.1 for k in G2_integers))

# Check 46: overshoot r-1 is close to 16*alpha/pi but not equal
r_minus_1 = r - 1
approx_16a_pi = 16 * float(alpha) / math.pi
check(46, "overshoot r-1 close to 16*alpha/pi (within 2%) but not equal",
      abs_err(r_minus_1, approx_16a_pi) / approx_16a_pi < 0.02
      and abs_err(r_minus_1, approx_16a_pi) > 1e-6)


# ===========================================================
# §7 — Verdict checks
# ===========================================================

# Check 47: OPEN verdict (no candidate achieves gap < 1e-5)
check(47, "OPEN verdict (no candidate achieves gap < 1e-5)", gap(c3) > 1e-5)

# Check 48: OPEN with candidate (best within 5% relative gap)
check(48, "OPEN with candidate (best within 5% relative gap)", gap(c3) / m_ratio < 0.05)

# Check 49: c3 relative gap ~3.67e-5
check(49, "c3 relative gap ~3.67e-5", rel_err(gap(c3)/m_ratio, 3.67e-5) < 0.05)

# Check 50: Full summary table — c3 strictly best among top 5
all_gaps = sorted([gap(c1), gap(c2), gap(c3),
                   gap(c_alpha_4pi), gap(c_cos30)])
check(50, "c3 strictly best among top 5 candidates", all_gaps[0] == gap(c3))


print(f"  ALPHA_INV = {float(ALPHA_INV):.9f}")
print(f"  alpha     = {float(alpha):.9e}")
print(f"  Delta     = {Delta:.9e}")
print(f"  alpha/2pi = {a_2pi:.9e}")
print(f"  r         = {r:.9f}")
print(f"  c* (best) = {c3:.9e}")
print(f"  gap(c*)   = {gap(c3):.6e}")
print(f"  rel gap   = {gap(c3)/m_ratio:.6e}")
print(f"  epsilon   = {epsilon:.6e}")
print(f"  N         = {N:.6f}")
print(f"  STATUS    : OPEN with candidate (best gap > 1e-5)")
print(f"\n{'='*60}\nRESULT: {PASS} PASS / {FAIL} FAIL")
import sys
sys.exit(0 if FAIL == 0 else 1)