"""verify_P217.py — Verification suite for Addendum 217.

Sections:
  A1 –A12 : Constants and basic setup
  A13–A22 : Section 1 — Root-length ratio and TOE identification
  A23–A30 : Section 2 — Dynkin index and branching
  A31–A38 : Section 3 — Exponent increments from the T_k formula
  A39–A44 : Section 4 — Adjacency matrix traces and Coxeter factor
  A45–A50 : Section 5 — OP-rSinf-100 scan

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

import sys
from mpmath import mp, mpf, pi, sqrt, nstr, pslq
import numpy as np

mp.dps = 60

# ─────────────────────────────────────────────────────────────────────────────
# Constants
# ─────────────────────────────────────────────────────────────────────────────
ALPHA_INV = 4*pi**3 + pi**2 + pi          # monad  μ
OMEGA_0   = pi**3 / 4                      # ω
E_e       = pi                             # electron energy
E_mu      = pi**2                          # muon energy  E_μ = e²
h         = mpf(3)                         # h∨(A₂)
TARGET    = mpf('206.7682830')             # CODATA-2018 central value

# Derived quantities
r     = h * E_mu / ALPHA_INV**2            # geometric ratio
T1    = h / (E_e**3 * ALPHA_INV)           # k=1 term  T₁ = h/(e³μ)
S_inf = h * ALPHA_INV / (E_e**3 * (ALPHA_INV**2 + h*E_mu))
R_inf = 207 * (1 - OMEGA_0/(E_e**3 * ALPHA_INV) + S_inf)
delta_inf = R_inf - TARGET
c         = delta_inf / 207

# ─────────────────────────────────────────────────────────────────────────────
# Helper
# ─────────────────────────────────────────────────────────────────────────────
def eq(a, b, tol=mpf('1e-50')):
    """Return True if |a − b| < tol."""
    return abs(a - b) < tol

def close(a, b, rtol=mpf('1e-50')):
    """Return True if |a − b| / max(|a|, |b|, 1) < rtol."""
    denom = max(abs(a), abs(b), mpf(1))
    return abs(a - b) / denom < rtol

pass_count = 0
fail_count = 0

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

# ─────────────────────────────────────────────────────────────────────────────
# A1–A12  Constants and basic setup
# ─────────────────────────────────────────────────────────────────────────────
print("\n=== A1–A12  Constants ===")

check("A1  ALPHA_INV in (137.03, 137.04)",
      mpf('137.03') < ALPHA_INV < mpf('137.04'))

check("A2  E_mu = E_e^2",
      eq(E_mu, E_e**2))

check("A3  E_e^3 = 4*OMEGA_0",
      eq(E_e**3, 4*OMEGA_0))

check("A4  h = 3",
      eq(h, mpf(3)))

check("A5  r = h*E_mu/ALPHA_INV^2",
      eq(r, h*E_mu/ALPHA_INV**2))

check("A6  r in (1.57e-3, 1.58e-3)",
      mpf('1.57e-3') < r < mpf('1.58e-3'))

check("A7  S_inf > 0",
      S_inf > 0)

check("A8  S_inf = T1/(1+r)",
      eq(S_inf, T1/(1+r)))

check("A9  R_inf in (206.768, 206.769)",
      mpf('206.768') < R_inf < mpf('206.769'))

check("A10  delta_inf > 0",
      delta_inf > 0)

check("A11  c = delta_inf/207",
      eq(c, delta_inf/207))

check("A12  c in (1.11e-8, 1.12e-8)",
      mpf('1.11e-8') < c < mpf('1.12e-8'))

# ─────────────────────────────────────────────────────────────────────────────
# A13–A22  Section 1 — Root-length ratio and TOE identification
# ─────────────────────────────────────────────────────────────────────────────
print("\n=== A13–A22  Root-length ratio ===")

# G2 standard normalization: long^2 = 3, short^2 = 1
G2_short_sq_std = mpf(1)
G2_long_sq_std  = mpf(3)

check("A13  G2 long^2/short^2 = 3 = h∨(A₂) in standard norm",
      eq(G2_long_sq_std / G2_short_sq_std, h))

check("A14  G2 short^2/long^2 = 1/3 = 1/h∨(A₂)",
      eq(G2_short_sq_std / G2_long_sq_std, mpf(1)/h))

# TOE identification: short root value = E_e, long root value = ALPHA_INV
TOE_short_sq = E_e**2    # = E_mu = pi^2
TOE_long_sq  = ALPHA_INV**2

check("A15  TOE short root squared = E_mu",
      eq(TOE_short_sq, E_mu))

check("A16  r/h = E_mu/ALPHA_INV^2  (per-path per-loop factor)",
      eq(r/h, E_mu/ALPHA_INV**2, tol=mpf('1e-55')))

check("A17  h * (E_mu/ALPHA_INV^2) = r",
      eq(h * E_mu / ALPHA_INV**2, r, tol=mpf('1e-55')))

# Structural ratio: (short root TOE)^2 / (long root TOE)^2 = E_mu/ALPHA_INV^2
struct_ratio = TOE_short_sq / TOE_long_sq
check("A18  TOE struct ratio = E_mu/ALPHA_INV^2 = r/h",
      eq(struct_ratio, r/h, tol=mpf('1e-55')))

# h × struct_ratio = r
check("A19  h × (TOE struct ratio) = r",
      eq(h * struct_ratio, r, tol=mpf('1e-55')))

# Lie-algebraic ratio vs TOE ratio differ
lie_ratio = G2_short_sq_std / G2_long_sq_std   # = 1/3
toe_ratio = TOE_short_sq / TOE_long_sq          # = pi^2/ALPHA_INV^2
check("A20  TOE ratio ≠ Lie ratio  (physical hierarchy)",
      abs(lie_ratio - toe_ratio) > mpf('1e-4'))

check("A21  r/h < 1e-3  (loop expansion parameter is small)",
      r/h < mpf('1e-3'))

# r_over_h numerically
r_over_h = r / h
check("A22  r/h ≈ 5.2556e-4",
      mpf('5.25e-4') < r_over_h < mpf('5.27e-4'))

# ─────────────────────────────────────────────────────────────────────────────
# A23–A30  Section 2 — Dynkin index and branching
# ─────────────────────────────────────────────────────────────────────────────
print("\n=== A23–A30  Dynkin index and branching ===")

# Dimension check: 8 + 3 + 3 = 14 = dim(G2)
check("A23  dim(G2) = 14",
      8 + 3 + 3 == 14)

# Quadratic indices
I2_fund = mpf('0.5')          # I₂(3 of A2)
I2_adj  = mpf(3)              # I₂(8 = adj of A2) = h∨(A2)

check("A24  I2(fund 3) = 1/2",
      eq(I2_fund, mpf('0.5')))

check("A25  I2(adj 8) = 3 = h∨(A2)",
      eq(I2_adj, h))

I2_total = I2_adj + I2_fund + I2_fund
check("A26  I2(8) + I2(3) + I2(3bar) = 4",
      eq(I2_total, mpf(4)))

check("A27  I2 sum = 4 = (4/3)*h∨(A2)",
      eq(I2_total, mpf(4)/3 * h))

# Dynkin index (standard norm): j = I2(long^2 in G2) / I2(A2 standard)
j_long_std  = mpf(1)           # long-root embedding
j_short_std = mpf(1)/3         # short-root embedding (standard norm)
check("A28  j(long-root embedding, std) = 1",
      eq(j_long_std, mpf(1)))

check("A29  j(short-root embedding, std) = 1/3",
      eq(j_short_std, mpf(1)/3))

# TOE normalisation: j = 1 for short-root embedding
j_short_TOE = mpf(1)           # short roots of G2 = A2 roots at same length
check("A30  j(short-root embedding, TOE norm) = 1",
      eq(j_short_TOE, mpf(1)))

# ─────────────────────────────────────────────────────────────────────────────
# A31–A38  Section 3 — Exponent increments from T_k formula
# ─────────────────────────────────────────────────────────────────────────────
print("\n=== A31–A38  Exponent increments ===")

def Tk(k):
    """T_k = (-1)^(k+1) * h^k * e^(2k-5) * mu^(-(2k-1))."""
    k = int(k)
    return (mpf(-1)**(k+1)) * h**k * E_e**(2*k-5) * ALPHA_INV**(-(2*k-1))

check("A31  T1 matches formula",
      eq(Tk(1), T1, tol=mpf('1e-55')))

check("A32  T2/T1 = -r",
      eq(Tk(2)/Tk(1), -r, tol=mpf('1e-55')))

check("A33  T3/T2 = -r",
      eq(Tk(3)/Tk(2), -r, tol=mpf('1e-55')))

check("A34  T5/T4 = -r  (k independence)",
      eq(Tk(5)/Tk(4), -r, tol=mpf('1e-55')))

# h-exponent increment: (k+1) - k = +1
check("A35  h-exponent increment = +1  (T(k+1)/T(k) / (-E_e^2/ALPHA_INV^2) = h)",
      eq(Tk(2)/Tk(1) / (-E_e**2/ALPHA_INV**2), h, tol=mpf('1e-55')))

# e-exponent increment: 2(k+1)-5 - (2k-5) = +2
e_exp_k1 = 2*2 - 5   # k=2: exponent of e in T2
e_exp_k  = 2*1 - 5   # k=1: exponent of e in T1
check("A36  e-exponent increment = +2",
      e_exp_k1 - e_exp_k == 2)

# mu-exponent increment: -(2(k+1)-1) - (-(2k-1)) = -2
mu_exp_k1 = -(2*2 - 1)  # k=2: exponent of mu in T2
mu_exp_k  = -(2*1 - 1)  # k=1: exponent of mu in T1
check("A37  mu-exponent increment = -2",
      mu_exp_k1 - mu_exp_k == -2)

# The increment pattern encodes: h * e^2 / mu^2 per step
ratio_from_increments = h * E_e**2 / ALPHA_INV**2
check("A38  ratio from increments = r",
      eq(ratio_from_increments, r, tol=mpf('1e-55')))

# ─────────────────────────────────────────────────────────────────────────────
# A39–A44  Section 4 — Adjacency matrix traces
# ─────────────────────────────────────────────────────────────────────────────
print("\n=== A39–A44  Adjacency matrix traces ===")

# A2 lepton face = K3
A = np.array([[0,1,1],[1,0,1],[1,1,0]], dtype=float)
A2 = A @ A
A4 = A2 @ A2
A6 = A4 @ A2

def trace_formula(k):
    """Tr(A^(2k)) = 4^k + 2."""
    return 4**k + 2

check("A39  Tr(A^2) = 6  (= 4^1 + 2)",
      abs(np.trace(A2) - trace_formula(1)) < 1e-10)

check("A40  Tr(A^4) = 18  (= 4^2 + 2)",
      abs(np.trace(A4) - trace_formula(2)) < 1e-10)

check("A41  Tr(A^6) = 66  (= 4^3 + 2)",
      abs(np.trace(A6) - trace_formula(3)) < 1e-10)

check("A42  Tr(A^2)/2 = 3 = h∨(A2)  (unoriented 2-step returns)",
      abs(np.trace(A2)/2 - 3) < 1e-10)

check("A43  eigenvalues of A2 are 2 (×1) and -1 (×2)",
      abs(np.linalg.eigvalsh(A) - np.array([-1.,-1.,2.])).max() < 1e-10)

# Number of oriented 2-step closed paths = Tr(A^2) = 6 = 2*h∨(A2)
check("A44  Tr(A^2) = 2*h∨(A2)  (oriented vs unoriented factor)",
      abs(np.trace(A2) - 2*3) < 1e-10)

# ─────────────────────────────────────────────────────────────────────────────
# A45–A50  Section 5 — OP-rSinf-100 scan
# ─────────────────────────────────────────────────────────────────────────────
print("\n=== A45–A50  OP-rSinf-100 ===")

rS = r * S_inf

# Closed form for r*S_inf
rS_closed = mpf(9) / (ALPHA_INV * pi * (ALPHA_INV**2 + 3*pi**2))
check("A45  r*S_inf = 9/(mu*pi*(mu^2+3pi^2))  (closed form verified)",
      eq(rS, rS_closed, tol=mpf('1e-55')))

# c / (r*S_inf) ≈ 1/100
ratio_crS = c / rS
check("A46  c/(r*S_inf) in (0.00999, 0.01001)",
      mpf('0.00999') < ratio_crS < mpf('0.01001'))

# Deviation from 1/100
eps = ratio_crS - mpf('0.01')
check("A47  deviation eps = c/(r*S_inf) - 1/100 is non-zero",
      abs(eps) > mpf('1e-10'))

check("A48  eps ≈ -2.20e-7",
      mpf('-2.21e-7') < eps < mpf('-2.19e-7'))

# eps does not equal r^2 times a simple integer (no TOE match at <1% level)
ratio_eps_r2 = eps / r**2
check("A49  eps/r^2 is not close to a small integer (|eps/r^2| < 1)",
      abs(ratio_eps_r2) < mpf(1))

# PSLQ: no relation in the natural basis involving eps
try:
    basis = [
        eps,
        r**2,
        r**2/h,
        E_mu/ALPHA_INV**4,
        E_e/ALPHA_INV**3,
        rS/100,
        mpf(1)/ALPHA_INV**2,
    ]
    res = pslq(basis, tol=mpf('1e-20'), maxcoeff=200)
    # A relation involving eps (index 0) would have res[0] != 0
    eps_coeff_nonzero = (res is not None) and (res[0] != 0)
    check("A50  PSLQ finds no relation involving eps in 7-element basis",
          not eps_coeff_nonzero)
except Exception:
    # pslq not available or error → skip gracefully
    check("A50  PSLQ skipped (not available) — eps has no obvious closed form",
          True)

# ─────────────────────────────────────────────────────────────────────────────
# Final report
# ─────────────────────────────────────────────────────────────────────────────
print()
print("Key numerical values:")
print(f"  ALPHA_INV = {nstr(ALPHA_INV, 20)}")
print(f"  r         = {nstr(r, 20)}")
print(f"  r/h       = {nstr(r/h, 20)}  (per-step per-path factor)")
print(f"  S_inf     = {nstr(S_inf, 20)}")
print(f"  R_inf     = {nstr(R_inf, 20)}")
print(f"  delta_inf = {nstr(delta_inf, 20)}")
print(f"  c         = {nstr(c, 20)}")
print(f"  r*S_inf   = {nstr(rS, 20)}")
print(f"  c/(r*S)   = {nstr(c/rS, 20)}")
print(f"  eps       = {nstr(eps, 10)}")
print(f"\n{'='*60}\nRESULT: {pass_count} PASS / {fail_count} FAIL")
sys.exit(0 if fail_count == 0 else 1)