"""
verify_P200.py — Verification suite for Addendum 200
G2 closed-form structure of the muon-to-electron mass ratio.

35+ assertions covering:
  - Constants and expansion coefficients
  - G2 structural identity (Proposition 2.1)
  - Tree-level decomposition of 207
  - Closed-form candidates from Table 1
  - Pade [1/1] construction and G2 expression of b1 (Proposition 5.1)
  - Four-loop residual from P199
  - Reorganisation identity: -N*(a/pi)^3 / a^2 = 11/16

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

from mpmath import mp, mpf, pi, sqrt, cos, exp, sin, besselj, fabs, nstr
mp.dps = 50

# ── Fundamental constants ────────────────────────────────────────────────────

ALPHA_INV = 4*pi**3 + pi**2 + pi
alpha     = 1 / ALPHA_INV
Omega0    = pi**3 / 4
R         = mpf('206.7682830')   # experimental m_mu/m_e
tree      = mpf('207')           # G2 tree-level

# ── G2 / A2 Lie data ─────────────────────────────────────────────────────────
dim_G2 = mpf('14')
hv_G2  = mpf('4')    # dual Coxeter number of G2
hv_A2  = mpf('3')    # dual Coxeter number of A2
dim_A2 = mpf('8')    # dimension of A2 (= su(3))

# ── Expansion coefficients ───────────────────────────────────────────────────
c1 = -1 / (2*pi)                  # coefficient of alpha (one-loop)
c2_pi = 1 / pi**2                 # "two-loop" piece of c2
c2_11 = mpf('11') / 16            # "three-loop reorganised" piece of c2
c2    = c2_pi + c2_11             # total O(alpha^2) coefficient
f     = c2                        # f from the paper

tol_tight  = mpf('1e-40')         # symbolic identity tolerance
tol_approx = mpf('1e-5')          # numerical approximation tolerance

PASS = FAIL = 0

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

# ============================================================
# §1 — Basic constants
# ============================================================

ok("alpha_inv > 137 and < 138",
   137 < ALPHA_INV < 138)

ok("alpha_inv close to 137.0363",
   fabs(ALPHA_INV - mpf('137.0363')) < mpf('0.0001'))

ok("alpha = 1/alpha_inv",
   fabs(alpha - 1/ALPHA_INV) < tol_tight)

ok("alpha ~ 7.297e-3",
   fabs(alpha - mpf('7.297e-3')) < mpf('1e-6'))

ok("Omega0 = pi^3/4",
   fabs(Omega0 - pi**3/4) < tol_tight)

ok("Omega0 ~ 7.752",
   fabs(Omega0 - mpf('7.752')) < mpf('0.001'))

ok("R = 206.7682830",
   fabs(R - mpf('206.7682830')) < tol_tight)

ok("tree = 207",
   fabs(tree - 207) < tol_tight)

# ============================================================
# §2 — Expansion coefficients
# ============================================================

ok("c1 = -1/(2pi)",
   fabs(c1 - (-1/(2*pi))) < tol_tight)

ok("c1 ~ -0.15915",
   fabs(c1 + mpf('0.15915')) < mpf('1e-4'))

ok("c2_pi = 1/pi^2",
   fabs(c2_pi - 1/pi**2) < tol_tight)

ok("c2_pi ~ 0.10132",
   fabs(c2_pi - mpf('0.10132')) < mpf('1e-4'))

ok("c2_11 = 11/16",
   fabs(c2_11 - mpf('11')/16) < tol_tight)

ok("c2_11 = 0.6875 exactly",
   fabs(c2_11 - mpf('0.6875')) < tol_tight)

ok("c2 = 1/pi^2 + 11/16",
   fabs(c2 - (1/pi**2 + mpf('11')/16)) < tol_tight)

ok("c2 ~ 0.78882",
   fabs(c2 - mpf('0.78882')) < mpf('1e-4'))

ok("c2 > 0",
   c2 > 0)

# ============================================================
# §3 — G2 structural identity (Proposition 2.1)
# ============================================================

f_G2 = 1/pi**2 + (dim_G2 - hv_A2) / hv_G2**2

ok("G2 identity: (dim_G2 - hv_A2)/hv_G2^2 = 11/16",
   fabs((dim_G2 - hv_A2) / hv_G2**2 - mpf('11')/16) < tol_tight)

ok("G2 identity: f = 1/pi^2 + (dim_G2-hv_A2)/hv_G2^2",
   fabs(f_G2 - f) < tol_tight)

ok("G2 identity matches c2",
   fabs(f_G2 - c2) < tol_tight)

ok("hv_G2^2 = 16",
   fabs(hv_G2**2 - 16) < tol_tight)

ok("dim_G2 - hv_A2 = 11",
   fabs(dim_G2 - hv_A2 - 11) < tol_tight)

# ============================================================
# §4 — Tree-level decomposition (§4 of paper)
# ============================================================

ok("207 = dim_G2*(dim_G2+1) - hv_A2",
   fabs(dim_G2*(dim_G2+1) - hv_A2 - 207) < tol_tight)

ok("207 = 14*15 - 3",
   14*15 - 3 == 207)

ok("207 = 9 * 23",
   9 * 23 == 207)

ok("23 = dim_G2 + dim_A2 + 1",
   fabs(dim_G2 + dim_A2 + 1 - 23) < tol_tight)

ok("hv_A2^2 = 9",
   fabs(hv_A2**2 - 9) < tol_tight)

# ============================================================
# §5 — Reorganisation identity
# ============================================================

N = -Omega0 * ALPHA_INV * (mpf('11')/4)

ok("N < 0",
   N < 0)

ok("N ~ -2921",
   fabs(N + 2921) < 1)

ok("-N*(a/pi)^3 / a^2 = 11/16  (reorganisation)",
   fabs((-N)*(alpha/pi)**3 / alpha**2 - mpf('11')/16) < tol_tight)

ok("-N*(a/pi)^3 / a^2 = c2_11",
   fabs((-N)*(alpha/pi)**3 / alpha**2 - c2_11) < tol_tight)

# ============================================================
# §6 — P199 formula and four-loop residual
# ============================================================

val_P199 = tree * (1 - alpha/(2*pi) + (alpha/pi)**2 - N*(alpha/pi)**3)
gap_P199 = fabs(val_P199 - R)

ok("P199 formula equals series to 45 digits",
   fabs(val_P199 - tree*(1 + c1*alpha + c2*alpha**2)) < mpf('1e-44'))

ok("P199 gap < 1e-5",
   gap_P199 < mpf('1e-5'))

ok("P199 gap ~ 8.81e-7",
   fabs(gap_P199 - mpf('8.81e-7')) < mpf('1e-8'))

ok("P199 formula > R  (series overshoots)",
   val_P199 > R)

# ============================================================
# §7 — Alpha-series (reorganised)
# ============================================================

val_series = tree * (1 + c1*alpha + c2*alpha**2)
gap_series = fabs(val_series - R)

ok("alpha-series = 206.76828388...",
   fabs(val_series - mpf('206.76828388')) < mpf('1e-6'))

ok("alpha-series gap < 1e-5",
   gap_series < tol_approx)

ok("alpha-series gap ~ 8.81e-7",
   fabs(gap_series - mpf('8.81e-7')) < mpf('1e-8'))

ok("alpha-series relative gap < 1e-8",
   gap_series/R < mpf('1e-8'))

# ============================================================
# §8 — Closed-form candidates (Table 1)
# ============================================================

# Candidate 2: exp(-a/2pi + (f - 1/(8pi^2))*a^2)
cand2 = exp(-alpha/(2*pi) + (f - 1/(8*pi**2))*alpha**2)
gap2 = fabs(tree*cand2 - R)
ok("candidate 2 gap < 1e-4",
   gap2 < mpf('1e-4'))
ok("candidate 2 gap ~ 8.9e-6",
   fabs(gap2 - mpf('8.9e-6')) < mpf('5e-7'))

# Candidate 8 (GON lapse): sqrt(1-alpha^2)
cand8 = sqrt(1 - alpha**2)
gap8 = fabs(tree*cand8 - R)
ok("candidate 8 (GON lapse) gap > 0.2",
   gap8 > mpf('0.2'))

# Best candidate is the polynomial — verify it beats all others
ok("polynomial beats candidate 2",
   gap_series < gap2)

ok("polynomial achieves sub-ppm relative accuracy",
   gap_series/R < mpf('1e-8'))

# ============================================================
# §9 — Padé [1/1] analysis
# ============================================================

b1 = 2*pi*c2
a1 = b1 - 1/(2*pi)

ok("b1 = 2*pi*c2",
   fabs(b1 - 2*pi*c2) < tol_tight)

ok("b1 ~ 4.9563",
   fabs(b1 - mpf('4.9563')) < mpf('1e-3'))

ok("a1 = b1 - 1/(2pi)",
   fabs(a1 - (b1 - 1/(2*pi))) < tol_tight)

ok("a1 ~ 4.7972",
   fabs(a1 - mpf('4.7972')) < mpf('1e-3'))

# G2 expression for b1 (Proposition 5.1)
b1_G2 = 2/pi + (dim_G2 - hv_A2)*pi/(2*hv_G2)
ok("b1 G2 form: 2/pi + (dim_G2-hv_A2)*pi/(2*hv_G2)",
   fabs(b1_G2 - b1) < tol_tight)

ok("(dim_G2-hv_A2)*pi/(2*hv_G2) = 11pi/8",
   fabs((dim_G2-hv_A2)*pi/(2*hv_G2) - 11*pi/8) < tol_tight)

# Evaluate Padé [1/1]
pade11 = (1 + a1*alpha)/(1 + b1*alpha)
val_pade = tree * pade11
gap_pade = fabs(val_pade - R)

ok("Pade[1/1] denominator > 1",
   1 + b1*alpha > 1)

ok("Pade[1/1] numerator > 1",
   1 + a1*alpha > 1)

ok("Pade[1/1] numerator < Pade[1/1] denominator",
   1 + a1*alpha < 1 + b1*alpha)

ok("Pade[1/1] gap ~ 3e-4",
   fabs(gap_pade - mpf('3.03e-4')) < mpf('5e-6'))

ok("Pade[1/1] gap > alpha-series gap (polynomial wins)",
   gap_pade > gap_series)

# ============================================================
# §10 — Verdict consistency
# ============================================================

ok("No candidate achieves gap < 1e-7  (verdict: OPEN)",
   gap_series > mpf('1e-7') and gap_pade > mpf('1e-7'))

ok("Best known gap is 8.81e-7 (four-loop residual)",
   fabs(gap_series - mpf('8.81e-7')) < mpf('1e-8'))

# ============================================================
# Summary
# ============================================================
print(f"\n{'='*55}")
print(f"  All {PASS} assertions PASSED")
print(f"  alpha        = {nstr(alpha, 12)}")
print(f"  c2           = {nstr(c2, 12)}")
print(f"  series value = {nstr(val_series, 12)}")
print(f"  series gap   = {nstr(gap_series, 4)}")
print(f"  Pade[1/1] gap= {nstr(gap_pade, 4)}")
print(f"  Verdict      : OPEN")
print(f"{'='*55}")

print(f"\n{'='*60}\nRESULT: {PASS} PASS / {FAIL} FAIL")