"""
verify_P204.py — Verification suite for Addendum P204
Loop Coefficient Self-Similarity and the Simplicial f-Vector.

All assertions use mpmath at dps=60. Run with:
    python verify_P204.py
All assertions must pass.

Sections:
  §1  TOE constants
  §2  Loop expansion coefficients
  §3  Ratios between consecutive coefficients
  §4  Self-similarity: |r34 + ALPHA_INV| / ALPHA_INV < 1e-3
  §5  Partial sums and CODATA straddling
  §6  Exact identity S4_hyp = S2 (Proposition 3.1)
  §7  f-vector of the regular 4-simplex via itertools.combinations
  §8  Binomial coefficient verification
  §9  E8 Coxeter number h(E8) = 30

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

import itertools
from math import comb

import mpmath
mpmath.mp.dps = 60

PASS = FAIL = 0
_N = 0


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


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

# ── CODATA 2018 target ────────────────────────────────────────────────────────
m_mu_me = mpmath.mpf('206.7682830')

# ── Loop coefficients ─────────────────────────────────────────────────────────
a1 = -1 / (2*pi)                                    # one-loop
a2 = 1/pi**2 + mpmath.mpf(11)/16                    # two-loop = f2
N  = -(mpmath.mpf(11)/4) * OMEGA_0 * ALPHA_INV      # three-loop N = a3
a3 = N

# Four-loop coefficient from CODATA inversion (P201)
f2  = a2
R3  = 207*(1 - alpha/(2*pi) + f2*alpha**2 + N*alpha**3)
c4  = (m_mu_me - R3) / (207 * alpha**4)
a4  = c4

# ── Ratios ────────────────────────────────────────────────────────────────────
r12 = a2 / a1
r23 = a3 / a2
r34 = a4 / a3

print("=" * 65)
print("verify_P204.py — P204 loop self-similarity verification")
print("=" * 65)

# ══════════════════════════════════════════════════════════════════════════════
# §1  TOE constants
# ══════════════════════════════════════════════════════════════════════════════

print("S1  TOE constants")

check(f"137 < ALPHA_INV = {float(ALPHA_INV)} < 138",
      137 < ALPHA_INV < 138)

check("OMEGA_0 definition",
      abs(OMEGA_0 - pi**3/4) < mpmath.mpf('1e-50'))

check("alpha = 1/ALPHA_INV",
      abs(alpha - 1/ALPHA_INV) < mpmath.mpf('1e-50'))

check("ALPHA_INV * alpha = 1 exactly",
      ALPHA_INV * alpha == 1 or abs(ALPHA_INV * alpha - 1) < mpmath.mpf('1e-55'))

print(f"  ALPHA_INV = {mpmath.nstr(ALPHA_INV, 15)}")
print(f"  OMEGA_0   = {mpmath.nstr(OMEGA_0, 12)}")
print(f"  alpha     = {mpmath.nstr(alpha, 12)}")

# ══════════════════════════════════════════════════════════════════════════════
# §2  Loop expansion coefficients
# ══════════════════════════════════════════════════════════════════════════════

print("S2  Loop expansion coefficients")

# a1
check("a1 definition",
      abs(a1 - (-1/(2*pi))) < mpmath.mpf('1e-50'))
check("a1 < 0",
      a1 < 0)
check("a1 ~ -0.15915",
      abs(a1 + mpmath.mpf('0.15915')) < mpmath.mpf('1e-4'))

# a2
check("a2 definition",
      abs(a2 - (1/pi**2 + mpmath.mpf(11)/16)) < mpmath.mpf('1e-50'))
check("a2 > 0",
      a2 > 0)
check("a2 ~ 0.78882",
      abs(a2 - mpmath.mpf('0.78882')) < mpmath.mpf('1e-4'))

# a3 = N
check("a3 = N < 0",
      a3 < 0)
check("a3 = N ~ -2921",
      abs(a3 + 2921) < 1)
check("a3 = N = -(11/4)*Omega0*ALPHA_INV",
      abs(a3 - (-(mpmath.mpf(11)/4)*OMEGA_0*ALPHA_INV)) < mpmath.mpf('1e-45'))

# a4 = c4
check("a4 = c4 > 0",
      a4 > 0)
check("a4 ~ 400305.875",
      abs(a4 - mpmath.mpf('400305.875')) < mpmath.mpf('0.01'))
check("floor(a4) = 400305",
      int(mpmath.floor(a4)) == 400305)

print(f"  a1 = {mpmath.nstr(a1, 10)}")
print(f"  a2 = {mpmath.nstr(a2, 10)}")
print(f"  a3 = N = {mpmath.nstr(a3, 12)}")
print(f"  a4 = c4 = {mpmath.nstr(a4, 15)}")

# ══════════════════════════════════════════════════════════════════════════════
# §3  Ratios between consecutive coefficients
# ══════════════════════════════════════════════════════════════════════════════

print("S3  Ratios between consecutive coefficients")

# r12 = a2/a1 = -2*pi*f2 ~ -4.953
check("r12 < 0 (sign change a1->a2)",
      r12 < 0)
check("r12 ~ -4.953",
      abs(r12 + mpmath.mpf('4.953')) < mpmath.mpf('0.01'))
check("r12 = -2*pi*a2 exactly",
      abs(r12 - (-2*pi*a2)) < mpmath.mpf('1e-50'))

# r23 = a3/a2 = N/f2 ~ -3703
check("r23 < 0 (sign change a2->a3)",
      r23 < 0)
check("|r23| > 3000 (N carries alpha^-1)",
      abs(r23) > 3000)
check("|r23| < 5000",
      abs(r23) < 5000)
check("r23 ~ -3703",
      abs(r23 + 3703) < 5)

# r34 = a4/a3 = c4/N ~ -137.04
check("r34 < 0 (sign change a3->a4)",
      r34 < 0)
check("|r34| > 130",
      abs(r34) > 130)
check("|r34| < 140",
      abs(r34) < 140)
check("r34 ~ -137.03",
      abs(r34 + mpmath.mpf('137.03')) < mpmath.mpf('0.1'))

print(f"  r12 = a2/a1 = {mpmath.nstr(r12, 10)}")
print(f"  r23 = a3/a2 = {mpmath.nstr(r23, 10)}")
print(f"  r34 = a4/a3 = {mpmath.nstr(r34, 10)}")
print(f"  -ALPHA_INV  = {mpmath.nstr(-ALPHA_INV, 10)}")

# ══════════════════════════════════════════════════════════════════════════════
# §4  Self-similarity: |r34 - (-ALPHA_INV)| / ALPHA_INV < 1e-3
# ══════════════════════════════════════════════════════════════════════════════

print("S4  Self-similarity: |r34 - (-ALPHA_INV)| / ALPHA_INV < 1e-3")

rel_dev = abs(r34 - (-ALPHA_INV)) / ALPHA_INV

# Main assertion: well within 0.1%
check(f"Self-similarity: rel_dev = {float(rel_dev):.2e} < 1e-3",
      rel_dev < mpmath.mpf('1e-3'))

# Tighter: the actual deviation is ~3.7e-6 (< 10 ppm)
check(f"Self-similarity tighter bound: rel_dev = {float(rel_dev):.2e} < 1e-4",
      rel_dev < mpmath.mpf('1e-4'))

# Direction: r34 < -ALPHA_INV (i.e. c4/N < -ALPHA_INV, equivalently c4 > -N*ALPHA_INV is FALSE;
# actually from P201: c4 < -N*ALPHA_INV means r34 = c4/N > c4/(-|N|) ... let's just check sign)
# From P201: c4 - (-N*ALPHA_INV) = -1.501 < 0, so c4 < -N*ALPHA_INV = |N|*ALPHA_INV
# Thus c4/|N| < ALPHA_INV, i.e. -c4/N < ALPHA_INV... wait N is negative:
# r34 = c4/N, N<0, c4>0 => r34 < 0.
# |r34| = c4/|N| < ALPHA_INV => r34 > -ALPHA_INV (r34 is less negative than -ALPHA_INV)
# Verify:
check("r34 > -ALPHA_INV (c4/|N| < ALPHA_INV, i.e. c4 < |N|*ALPHA_INV, consistent with P201)",
      r34 > -ALPHA_INV)

# The deviation from near-coincidence (P201 eq. 7.4): |c4 + N*ALPHA_INV| < 2
nc_residual = c4 + N * ALPHA_INV          # = c4 - |N|*ALPHA_INV = c4 - (-N*ALPHA_INV)
check("c4 < -N*ALPHA_INV (P201 direction)",
      nc_residual < 0)
check("|c4 + N*ALPHA_INV| < 2 (P201 §7.4)",
      abs(nc_residual) < 2)
check("nc residual ~ -1.501 (P201)",
      abs(nc_residual + mpmath.mpf('1.501')) < mpmath.mpf('0.01'))

dev_pct = float(rel_dev) * 100
print(f"  rel_dev = {float(rel_dev):.4e}  ({dev_pct:.5f}%)")
print(f"  nc_residual (c4 + N*ALPHA_INV) = {mpmath.nstr(nc_residual, 10)}")

# ══════════════════════════════════════════════════════════════════════════════
# §5  Partial sums and CODATA straddling
# ══════════════════════════════════════════════════════════════════════════════

print("S5  Partial sums and CODATA straddling")

S2 = 207*(1 + a1*alpha + a2*alpha**2)
S3 = 207*(1 + a1*alpha + a2*alpha**2 + a3*alpha**3)
S4_exact = 207*(1 + a1*alpha + a2*alpha**2 + a3*alpha**3 + a4*alpha**4)

# Hypothesis value
a4_hyp = (-ALPHA_INV) * a3          # = -alpha^-1 * N = |N| * ALPHA_INV > 0
S4_hyp = 207*(1 + a1*alpha + a2*alpha**2 + a3*alpha**3 + a4_hyp*alpha**4)

# S2 overshoots CODATA (P200 result)
check("S2 > CODATA (two-loop overshoots)",
      S2 > m_mu_me)
check("S2 ~ 206.76828388",
      abs(S2 - mpmath.mpf('206.76828388')) < mpmath.mpf('1e-7'))
check("S2 - CODATA ~ 8.8e-7",
      abs(S2 - m_mu_me - mpmath.mpf('8.8e-7')) < mpmath.mpf('5e-9'))

# S3 undershoots CODATA (P201 result)
check("S3 < CODATA (three-loop undershoots)",
      S3 < m_mu_me)
check("S3 ~ 206.533",
      abs(S3 - mpmath.mpf('206.533')) < mpmath.mpf('0.001'))

# S4_exact = CODATA by construction
check("S4_exact = CODATA by construction",
      abs(S4_exact - m_mu_me) < mpmath.mpf('1e-9'))

# a4_hyp > 0 and > c4 (since near_c4 = -N*ALPHA_INV > c4)
check("a4_hyp > 0",
      a4_hyp > 0)
check("a4_hyp > c4 (hypothesis overshoots)",
      a4_hyp > c4)

# Straddling: S3 < CODATA < S4_hyp
check("S3 < CODATA (lower bound of straddle)",
      S3 < m_mu_me)
check("S4_hyp > CODATA (upper bound of straddle)",
      S4_hyp > m_mu_me)

print(f"  S2        = {mpmath.nstr(S2, 14)}")
print(f"  S3        = {mpmath.nstr(S3, 12)}")
print(f"  S4_exact  = {mpmath.nstr(S4_exact, 14)}")
print(f"  S4_hyp    = {mpmath.nstr(S4_hyp, 14)}")
print(f"  CODATA    = {mpmath.nstr(m_mu_me, 14)}")
print(f"  Straddle: S3 ({float(S3):.6f}) < CODATA < S4_hyp ({float(S4_hyp):.8f}): "
      f"{bool(S3 < m_mu_me < S4_hyp)}")

# ══════════════════════════════════════════════════════════════════════════════
# §6  Exact identity S4_hyp = S2 (Proposition 3.1)
# ══════════════════════════════════════════════════════════════════════════════

print("S6  Exact identity S4_hyp = S2 (Proposition 3.1)")

# Identity: S4_hyp = S2 exactly because ALPHA_INV * alpha = 1
# Proof (algebraic): S4_hyp = S3 + 207*a4_hyp*alpha^4
#                           = (S2 + 207*N*alpha^3) + 207*(-N*ALPHA_INV)*alpha^4
#                           = S2 + 207*N*alpha^3*(1 - ALPHA_INV*alpha)
#                           = S2 + 207*N*alpha^3*(1 - 1)
#                           = S2

diff_S4hyp_S2 = S4_hyp - S2
check(f"Proposition 3.1: S4_hyp = S2 exactly; diff = {float(diff_S4hyp_S2):.2e}",
      abs(diff_S4hyp_S2) < mpmath.mpf('1e-50'))

# Corollary: S4_hyp - CODATA = S2 - CODATA = delta2
delta2 = S2 - m_mu_me
check("S4_hyp - CODATA = S2 - CODATA",
      abs(S4_hyp - m_mu_me - delta2) < mpmath.mpf('1e-50'))

# delta2 ~ 8.8e-7 (P200 residual)
check("delta2 ~ 8.8e-7",
      abs(delta2 - mpmath.mpf('8.8e-7')) < mpmath.mpf('5e-9'))
check("delta2 > 0 (S2 overshoots)",
      delta2 > 0)

# The hypothesis overshoots CODATA by the same amount as the two-loop truncation
gap_hyp = S4_hyp - m_mu_me
check("hypothesis overshoot = two-loop overshoot (exact identity)",
      abs(gap_hyp - delta2) < mpmath.mpf('1e-50'))

print(f"  S4_hyp - S2     = {float(diff_S4hyp_S2):.2e}  (should be 0 to 1e-50)")
print(f"  S2 - CODATA     = {mpmath.nstr(delta2, 8)}")
print(f"  S4_hyp - CODATA = {mpmath.nstr(gap_hyp, 8)}")

# ══════════════════════════════════════════════════════════════════════════════
# §7  f-vector of the regular 4-simplex via itertools.combinations
# ══════════════════════════════════════════════════════════════════════════════

print("S7  f-vector of the regular 4-simplex via itertools.combinations")

VERTS = [0, 1, 2, 3, 4]                           # 5 vertices of the 4-simplex

faces_0 = list(itertools.combinations(VERTS, 1))   # vertices  (0-faces)
faces_1 = list(itertools.combinations(VERTS, 2))   # edges     (1-faces)
faces_2 = list(itertools.combinations(VERTS, 3))   # triangles (2-faces)
faces_3 = list(itertools.combinations(VERTS, 4))   # tetrahedra (3-faces)

check(f"5 vertices, got {len(faces_0)}",
      len(faces_0) == 5)
check(f"10 edges, got {len(faces_1)}",
      len(faces_1) == 10)
check(f"10 triangles, got {len(faces_2)}",
      len(faces_2) == 10)
check(f"5 tetrahedra, got {len(faces_3)}",
      len(faces_3) == 5)

f_vec = (len(faces_0), len(faces_1), len(faces_2), len(faces_3))
check(f"f-vector should be (5,10,10,5), got {f_vec}",
      f_vec == (5, 10, 10, 5))

total_faces = sum(f_vec)
check(f"total proper faces should be 30, got {total_faces}",
      total_faces == 30)

# Palindromic symmetry
check("f0 = f3 (palindrome)",
      f_vec[0] == f_vec[3])
check("f1 = f2 (palindrome)",
      f_vec[1] == f_vec[2])

print(f"  Vertices (0-faces): {len(faces_0)}")
print(f"  Edges    (1-faces): {len(faces_1)}")
print(f"  Triangles(2-faces): {len(faces_2)}")
print(f"  Tetrahedra(3-faces):{len(faces_3)}")
print(f"  f-vector = {f_vec}")
print(f"  Total proper faces = {total_faces}")

# ══════════════════════════════════════════════════════════════════════════════
# §8  Binomial coefficient verification
# ══════════════════════════════════════════════════════════════════════════════

print("S8  Binomial coefficient verification")

# k-faces of the 4-simplex (on 5 vertices): count = C(5, k+1)
check("C(5,1) = 5  (vertices)",
      comb(5, 1) == 5)
check("C(5,2) = 10 (edges)",
      comb(5, 2) == 10)
check("C(5,3) = 10 (triangles)",
      comb(5, 3) == 10)
check("C(5,4) = 5  (tetrahedra)",
      comb(5, 4) == 5)

# f-vector matches binomial coefficients
binom_fvec = tuple(comb(5, k) for k in range(1, 5))
check(f"f-vector should equal (C(5,1),..,C(5,4)), got {binom_fvec}",
      binom_fvec == (5, 10, 10, 5))

# Total = 30
total_binom = sum(comb(5, k) for k in range(1, 5))
check(f"sum C(5,k) k=1..4 = 30, got {total_binom}",
      total_binom == 30)

# Palindromic symmetry from binomial symmetry C(n,k) = C(n,n-k)
check("C(5,1) = C(5,4) (palindrome)",
      comb(5, 1) == comb(5, 4))
check("C(5,2) = C(5,3) (palindrome)",
      comb(5, 2) == comb(5, 3))

# Row-sum identity: sum_{k=0}^{5} C(5,k) = 2^5 = 32
# Proper faces exclude the empty face C(5,0)=1 and full simplex C(5,5)=1
check("2^5 = 32",
      sum(comb(5, k) for k in range(6)) == 32)
check("32 - 1 - 1 = 30 proper faces",
      32 - comb(5, 0) - comb(5, 5) == 30)

# The 5 Wheel operators label the vertices
N_WHEEL_OPS = 5  # Fork, Weld, Plateau, Oscillate, Perturb
check("5 Wheel operators (CLAUDE.md)",
      N_WHEEL_OPS == 5)
check("Wheel operators = C(5,1) = 5 vertices",
      N_WHEEL_OPS == comb(5, 1))

print(f"  (C(5,1), C(5,2), C(5,3), C(5,4)) = {binom_fvec}")
print(f"  sum = {total_binom}")

# ══════════════════════════════════════════════════════════════════════════════
# §9  E8 Coxeter number h(E8) = 30
# ══════════════════════════════════════════════════════════════════════════════

print("S9  E8 Coxeter number h(E8) = 30")

H_E8 = 30   # Coxeter number of E8 (classical; rank 8, dim 248, exponents 1,7,11,13,17,19,23,29)

check("h(E8) = 30",
      H_E8 == 30)
check("total f-vector = h(E8)",
      total_faces == H_E8)
check("binomial sum = h(E8)",
      total_binom == H_E8)

# Additional E8 invariants (constants, not derived)
DIM_E8  = 248
RANK_E8 = 8
check("dim E8 = 248",
      DIM_E8  == 248)
check("rank E8 = 8",
      RANK_E8 == 8)
# Coxeter number identity: sum of exponents + rank = rank * h
# Exponents of E8: 1,7,11,13,17,19,23,29. Sum = 120. rank*h - rank = 8*30 - 8 = 232. Hmm.
# Alternatively: number of positive roots = rank*(h)/2 = 8*30/2 = 120. Checked.
N_POS_ROOTS_E8 = 120
check("rank(E8) * h(E8) / 2 = 120 positive roots",
      RANK_E8 * H_E8 // 2 == N_POS_ROOTS_E8)

# Consistency: 30 = 5 + 10 + 10 + 5 (direct check)
check("5+10+10+5 = h(E8)",
      5 + 10 + 10 + 5 == H_E8)

print(f"  h(E8) = {H_E8}")
print(f"  f-vector sum = {total_faces} = h(E8) = {H_E8}: {total_faces == H_E8}")

# ══════════════════════════════════════════════════════════════════════════════
# Final summary
# ══════════════════════════════════════════════════════════════════════════════

print()
print("=" * 65)
print("  TOE constants:")
print(f"    ALPHA_INV = {mpmath.nstr(ALPHA_INV, 15)}")
print(f"    OMEGA_0   = {mpmath.nstr(OMEGA_0, 12)}")
print()
print("  Loop coefficients:")
print(f"    a1 = {mpmath.nstr(a1, 10)}")
print(f"    a2 = {mpmath.nstr(a2, 10)}")
print(f"    a3 = N = {mpmath.nstr(a3, 12)}")
print(f"    a4 = c4 = {mpmath.nstr(a4, 15)}")
print()
print("  Ratios:")
print(f"    r12 = a2/a1 = {mpmath.nstr(r12, 10)}")
print(f"    r23 = a3/a2 = {mpmath.nstr(r23, 10)}")
print(f"    r34 = a4/a3 = {mpmath.nstr(r34, 10)}")
print(f"    -ALPHA_INV  = {mpmath.nstr(-ALPHA_INV, 10)}")
print(f"    self-similarity deviation: {float(rel_dev*100):.6f}%  ({float(rel_dev*1e6):.2f} ppm)")
print()
print("  Partial sums:")
print(f"    S2        = {mpmath.nstr(S2, 14)}")
print(f"    S3        = {mpmath.nstr(S3, 12)}")
print(f"    CODATA    = {mpmath.nstr(m_mu_me, 14)}")
print(f"    S4_exact  = {mpmath.nstr(S4_exact, 14)}")
print(f"    S4_hyp    = {mpmath.nstr(S4_hyp, 14)}  (= S2 exactly)")
print(f"    Straddle: S3 < CODATA < S4_hyp = S2: {bool(S3 < m_mu_me < S4_hyp)}")
print()
print("  Simplicial f-vector (4-simplex):")
print(f"    f-vector = {f_vec}  (= (C(5,1),..,C(5,4)))")
print(f"    total proper faces = {total_faces} = h(E8) = {H_E8}")
print()
print("  Verdict: OPEN — self-similarity established numerically (3.7e-6 rel.)")
print("           Propositions 3.1 and 5.1 exact.")
print("           OP-C4 (P201) remains open.")
print("=" * 65)

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