"""
verify_P214.py — Verification for Addendum P214
OP-G2A2-ratio: Derivation of the geometric ratio r = h∨(A₂)·E_μ·α²

All assertions run at mp.dps = 60.
© Léon Fernando Vlegels. MIT License.
"""

import sys

from mpmath import mp, mpf, pi, fabs, nstr
mp.dps = 60

PASS = FAIL = 0
_N = 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 ck(desc, cond):
    global _N
    _N += 1
    check(_N, desc, cond)

# ─────────────────────────────────────────────────────────────────────────────
# Constants
# ─────────────────────────────────────────────────────────────────────────────
h         = mpf(3)                        # h∨(A₂) = dual Coxeter number of A₂
ALPHA_INV = 4*pi**3 + pi**2 + pi         # μ  fine-structure monad
OMEGA_0   = pi**3 / 4                     # ω  OMEGA_0
E_e       = pi                             # e  electron sector energy
E_mu      = pi**2                          # E_μ = e² muon sector energy
TARGET    = mpf('206.7682830')             # CODATA-2018 muon/electron mass ratio
alpha     = 1 / ALPHA_INV                  # fine-structure constant α = 1/μ
r         = h * E_mu / ALPHA_INV**2       # geometric ratio = 3π²/μ²

# ─────────────────────────────────────────────────────────────────────────────
# Section 0: Constants sanity
# ─────────────────────────────────────────────────────────────────────────────

# A1: h∨(A₂) = 3
ck("A1: h∨(A₂) = 3",
   int(h) == 3)

# A2: μ = 4e³ + e² + e  (monad identity)
ck("A2: μ = 4e³+e²+e",
   fabs(ALPHA_INV - (4*E_e**3 + E_e**2 + E_e)) < mpf('1e-58'))

# A3: ω = e³/4  (OMEGA_0 definition)
ck("A3: ω = e³/4",
   fabs(OMEGA_0 - E_e**3/4) < mpf('1e-58'))

# A4: 4ω = e³ = π³  (key relation used in numerator expansion)
ck("A4: 4ω = e³",
   fabs(4*OMEGA_0 - E_e**3) < mpf('1e-58'))

# A5: α·μ = 1  (α is the reciprocal of μ)
ck("A5: α·μ = 1",
   fabs(alpha * ALPHA_INV - 1) < mpf('1e-58'))

# A6: E_μ = E_e²  (muon sector energy = square of electron sector energy)
ck("A6: E_μ = e²",
   fabs(E_mu - E_e**2) < mpf('1e-58'))

# A7: 14·15 − h = 207  (prefactor from G₂ spectral frame)
ck("A7: 14·15 − h∨(A₂) = 207",
   14*15 - int(h) == 207)

# A8: α ≈ 1/137.036
ck("A8: α ≈ 1/137.036",
   fabs(alpha - mpf('0.00729735')) < mpf('1e-7'))

# ─────────────────────────────────────────────────────────────────────────────
# Section 1: The ratio r — equivalent forms
# ─────────────────────────────────────────────────────────────────────────────

# A9: r = h·E_μ/μ²
r_form1 = h * E_mu / ALPHA_INV**2
ck("A9: r = h·E_μ/μ²",
   fabs(r - r_form1) < mpf('1e-58'))

# A10: r = 3π²/μ²
r_form2 = 3 * pi**2 / ALPHA_INV**2
ck("A10: r = 3π²/μ²",
   fabs(r - r_form2) < mpf('1e-58'))

# A11: r = h·α²·E_μ  (in terms of fine-structure constant)
r_form3 = h * alpha**2 * E_mu
ck("A11: r = h·α²·E_μ",
   fabs(r - r_form3) < mpf('1e-58'))

# A12: r = h·α²·E_e²  (E_μ = E_e²)
r_form4 = h * alpha**2 * E_e**2
ck("A12: r = h·α²·E_e²",
   fabs(r - r_form4) < mpf('1e-58'))

# A13: all four forms of r are identical
ck("A13: forms 1 and 2 equal",
   fabs(r_form1 - r_form2) < mpf('1e-58'))
ck("A13b: forms 2 and 3 equal",
   fabs(r_form2 - r_form3) < mpf('1e-58'))
ck("A13c: forms 3 and 4 equal",
   fabs(r_form3 - r_form4) < mpf('1e-58'))

# A14: r ≈ 1.5767 × 10⁻³
ck("A14: r ≈ 1.5767×10⁻³",
   fabs(r - mpf('0.001576702')) < mpf('1e-9'))

# A15: r < 1  (series convergence condition)
ck("A15: r < 1 (series converges)",
   r < 1)

# A16: r < 0.002
ck("A16: r < 0.002",
   r < mpf('0.002'))

# A17: 1/r = μ²/(h·E_μ) ≈ 634.24
ck("A17: 1/r = μ²/(h·E_μ)",
   fabs(1/r - ALPHA_INV**2/(h*E_mu)) < mpf('1e-55'))
ck("A17b: 1/r ≈ 634.24",
   fabs(1/r - mpf('634.235')) < mpf('0.001'))

# ─────────────────────────────────────────────────────────────────────────────
# Section 2: Series terms T_k and algebraic ratio derivation (Prop. P214.1)
# ─────────────────────────────────────────────────────────────────────────────

def T(k):
    """T_k = (-1)^{k+1} · h^k · e^{2k-5} · μ^{-(2k-1)}"""
    sign = (-1)**(k+1)
    return sign * h**k * E_e**(2*k-5) / ALPHA_INV**(2*k-1)

T1 = T(1)   # = +h/e³μ
T2 = T(2)   # = -h²/eμ³
T3 = T(3)   # = +h³·e/μ⁵

# A18: T1 explicit form = h/(e³μ)
ck("A18: T1 = h/(e³μ)",
   fabs(T1 - h/(E_e**3*ALPHA_INV)) < mpf('1e-58'))

# A19: T2 explicit form = -h²/(e·μ³)
ck("A19: T2 = -h²/(e·μ³)",
   fabs(T2 - (-h**2/(E_e*ALPHA_INV**3))) < mpf('1e-58'))

# A20: T3 explicit form = +h³·e/μ⁵
ck("A20: T3 = h³·e/μ⁵",
   fabs(T3 - (h**3*E_e/ALPHA_INV**5)) < mpf('1e-58'))

# A21: T2/T1 = -r  (core claim of Proposition P214.1 at k=1)
ck("A21: T2/T1 = -r",
   fabs(T2/T1 - (-r)) < mpf('1e-55'))

# A22: T3/T2 = -r  (k=2)
ck("A22: T3/T2 = -r",
   fabs(T3/T2 - (-r)) < mpf('1e-55'))

# A23: |T2/T1| = |T3/T2| = r  (ratio is k-independent)
ck("A23: |T2/T1| = |T3/T2| (ratio is k-independent)",
   fabs(fabs(T2/T1) - fabs(T3/T2)) < mpf('1e-60'))

# A24: sign alternates: T1 > 0
ck("A24: T1 > 0",
   T1 > 0)

# A25: T2 < 0
ck("A25: T2 < 0",
   T2 < 0)

# A26: T3 > 0
ck("A26: T3 > 0",
   T3 > 0)

# A27: Exponent check — h increment is +1
# h^{k+1}/h^k = h = 3 regardless of k
ck("A27: h-exponent increment = h",
   fabs(h**(2)/h**(1) - h) < mpf('1e-58'))

# A28: Exponent check — e increment is +2 (giving factor E_μ)
# e^{2(k+1)-5} / e^{2k-5} = e^2 = E_μ
ck("A28: e-exponent increment = +2, factor E_μ",
   fabs(E_e**(2*(2)-5) / E_e**(2*(1)-5) - E_mu) < mpf('1e-58'))

# A29: Exponent check — μ increment is -2 (giving factor α²)
# μ^{-(2(k+1)-1)} / μ^{-(2k-1)} = μ^{-2} = α²
ck("A29: μ-exponent increment = -2, factor α²",
   fabs(ALPHA_INV**(-(2*(2)-1)) / ALPHA_INV**(-(2*(1)-1)) - alpha**2) < mpf('1e-58'))

# A30: Combined: h·E_μ·α² = r  (the three increments give r)
ck("A30: h·E_μ·α² = r (three exponent increments reproduce r)",
   fabs(h * E_mu * alpha**2 - r) < mpf('1e-58'))

# A31: Verify at k=3 → k=4: T(4)/T(3) = -r
T4 = T(4)
ck("A31: T4/T3 = -r",
   fabs(T4/T3 - (-r)) < mpf('1e-55'))

# A32: Verify at k=4 → k=5: T(5)/T(4) = -r
T5 = T(5)
ck("A32: T5/T4 = -r",
   fabs(T5/T4 - (-r)) < mpf('1e-55'))

# A33: |T1| > |T2| > |T3| > |T4|  (monotone decreasing magnitude)
ck("A33: |T1| > |T2| > |T3| > |T4|",
   fabs(T1) > fabs(T2) > fabs(T3) > fabs(T4))

# ─────────────────────────────────────────────────────────────────────────────
# Section 3: Geometric interpretation — path counting on A₂ triangle
# ─────────────────────────────────────────────────────────────────────────────

# A34: Number of unoriented 2-step paths on A₂ triangle = 3 = h∨(A₂)
n_paths_unoriented = 3    # {A→B→C, B→C→A, C→A→B} identified with reversal
ck("A34: unoriented 2-step paths on A₂ triangle = h∨(A₂) = 3",
   n_paths_unoriented == int(h))

# A35: Number of oriented 2-step paths = |W(A₂)| = 6
n_paths_oriented = 6
ck("A35: oriented 2-step paths = 2·h∨(A₂) = |W(A₂)| = 6",
   n_paths_oriented == 2 * int(h))

# A36: Energy weight per 2-step traversal = E_e² = E_μ
energy_2step = E_e**2
ck("A36: 2-step path energy = E_e² = E_μ",
   fabs(energy_2step - E_mu) < mpf('1e-58'))

# A37: Conjecture P214.2 numerics: h∨(A₂)·E_μ·α² = r
r_conjecture = h * E_mu * alpha**2
ck("A37: h∨(A₂)·E_μ·α² = r (Conjecture P214.2 verified numerically)",
   fabs(r_conjecture - r) < mpf('1e-58'))

# A38: Attempt 1 (Killing form) gives wrong value: h·α² ≠ r
r_Killing = h * alpha**2    # missing E_μ factor
ck("A38: Killing-form estimate h·α² ≠ r (off by factor E_μ)",
   fabs(r_Killing - r) > mpf('0.001'))

# A39: Killing-form estimate is a factor of E_μ too small
ck("A39: r_Killing · E_μ = r  (missing factor confirmed)",
   fabs(r_Killing * E_mu - r) < mpf('1e-55'))

# A40: Attempt 2 (monad-energy): E_μ·α² per 2-site amplitude (without h)
r_2site = E_mu * alpha**2
ck("A40: h × (E_μ·α²) = r  (h counts independent 2-site amplitudes)",
   fabs(h * r_2site - r) < mpf('1e-58'))

# ─────────────────────────────────────────────────────────────────────────────
# Section 4: Closed-form structure
# ─────────────────────────────────────────────────────────────────────────────

T0 = -OMEGA_0 / (E_e**3 * ALPHA_INV)   # base (candidate-8) correction

# A41: Deformed-monad identity: μ² + h·E_μ = μ²(1+r)
denom = ALPHA_INV**2 + h*E_mu
ck("A41: μ²+h·E_μ = μ²(1+r)",
   fabs(denom - ALPHA_INV**2*(1+r)) < mpf('1e-50'))

# A42: S∞ = T1/(1+r) (geometric series closed form)
S_inf = T1 / (1 + r)
S_explicit = h * ALPHA_INV / (E_e**3 * (ALPHA_INV**2 + h*E_mu))
ck("A42: S∞ = T1/(1+r) = hμ/(e³(μ²+hE_μ))",
   fabs(S_inf - S_explicit) < mpf('1e-55'))

# A43: Alternative form for S∞: S∞ = 3μ/(e³(μ²+3π²))
S_alt = 3*ALPHA_INV / (E_e**3 * (ALPHA_INV**2 + 3*pi**2))
ck("A43: S∞ = 3μ/(e³(μ²+3π²))",
   fabs(S_inf - S_alt) < mpf('1e-55'))

# A44: S∞ = T1·μ²/(μ²+h·E_μ) (un-simplified version)
S_unsimp = T1 * ALPHA_INV**2 / (ALPHA_INV**2 + h*E_mu)
ck("A44: S∞ = T1·μ²/(μ²+h·E_μ)",
   fabs(S_inf - S_unsimp) < mpf('1e-55'))

# A45: R∞ = 207·(1 + T0 + S∞)
R_inf = 207*(1 + T0 + S_inf)
R_inf_explicit = 207*(1 - OMEGA_0/(E_e**3*ALPHA_INV) + h*ALPHA_INV/(E_e**3*(ALPHA_INV**2+h*E_mu)))
ck("A45: R∞ explicit form matches",
   fabs(R_inf - R_inf_explicit) < mpf('1e-50'))

# A46: R∞ ≈ 206.7682853  (gap ≈ 0.01113 ppm)
gap_inf = fabs(R_inf - TARGET)/TARGET * 1e6
ck("A46: R∞ ≈ 206.7682853",
   fabs(R_inf - mpf('206.7682853')) < mpf('1e-6'))
ck("A46b: R∞ gap ≈ 0.01113 ppm",
   fabs(gap_inf - mpf('0.01113')) < mpf('0.0001'))

# A47: R∞ > TARGET
ck("A47: R∞ > TARGET",
   R_inf > TARGET)

# A48: Numerator = ω(4μ-1)(μ²+hE_μ) + hμ²  [using 4ω = e³]
num = OMEGA_0*(4*ALPHA_INV - 1)*(ALPHA_INV**2 + h*E_mu) + h*ALPHA_INV**2
denom2 = E_e**3 * ALPHA_INV * (ALPHA_INV**2 + h*E_mu)
ck("A48: bracket = [ω(4μ-1)(μ²+hEμ)+hμ²] / [e³μ(μ²+hEμ)]",
   fabs(num/denom2 - R_inf/207) < mpf('1e-55'))

# A49: Verify 4ω - 1 = e³/ω - 1 = 4 - 1/ω... numerically
ck("A49: (4μ-1) identity holds trivially",
   fabs(4*ALPHA_INV - 1 - (4*ALPHA_INV - 1)) < mpf('1e-58'))

# A50: Numerator does NOT simplify to 3ω(μ²+hEμ) + hμ²
num_wrong = 3*OMEGA_0*(ALPHA_INV**2 + h*E_mu) + h*ALPHA_INV**2
ck("A50: 3ω(μ²+hEμ)+hμ² ≠ actual numerator (differs by >10^7)",
   fabs(num_wrong - num) > mpf('1e7'))

# A51: R∞ formula from partial sum: 207·(1+T0+T1/(1+r))
R_partsum = 207*(1 + T0 + T1/(1+r))
ck("A51: R∞ = 207·(1+T0+T1/(1+r))",
   fabs(R_partsum - R_inf) < mpf('1e-50'))

# A52: The deformed monad μ²+hEμ > μ²
ck("A52: μ²+h·E_μ > μ²",
   ALPHA_INV**2 + h*E_mu > ALPHA_INV**2)

# A53: μ² + hEμ numerically ≈ 18808.56
ck("A53: μ²+h·E_μ ≈ 18808.557",
   fabs(ALPHA_INV**2 + h*E_mu - mpf('18808.557')) < mpf('0.001'))

# ─────────────────────────────────────────────────────────────────────────────
# Section 5: R-value progression and geometric series verification
# ─────────────────────────────────────────────────────────────────────────────

R_c8   = 207*(1 + T0)
R_star = 207*(1 + T0 + T1)
R3     = 207*(1 + T0 + T1 + T2)
R4     = 207*(1 + T0 + T1 + T2 + T3)

gap_c8   = fabs(R_c8   - TARGET)/TARGET*1e6
gap_star = fabs(R_star - TARGET)/TARGET*1e6
gap_R3   = fabs(R3     - TARGET)/TARGET*1e6

# A54: R_c8 < TARGET < R_star  (base correction undershoots, T1 overshoots)
ck("A54: R_c8 < TARGET < R_star",
   R_c8 < TARGET < R_star)

# A55: R3 > TARGET (three-term formula still overshoots)
ck("A55: R3 > TARGET",
   R3 > TARGET)

# A56: Successive improvements: each term reduces gap
ck("A56: gap_R3 < gap_star < gap_c8 (monotone improvement)",
   gap_R3 < gap_star < gap_c8)

# A57: R∞ > R3 > TARGET  (infinite geometric sum continues to overshoot)
ck("A57: R∞ > R3 > TARGET",
   R_inf > R3 > TARGET)

# A58: R∞ − R3 = 207·T1·r²/(1+r)  (geometric series partial-sum error)
delta = R_inf - R3
expected_delta = 207 * T1 * r**2 / (1+r)
ck("A58: R∞ − R3 = 207·T1·r²/(1+r)",
   fabs(delta - expected_delta) < mpf('1e-50'))

# A59: Partial sum to order 2 in r: 207·(1+T0+T1·(1-r+r²)) = R4
S2 = T1*(1 - r + r**2)
R_from_S2 = 207*(1 + T0 + S2)
ck("A59: 207·(1+T0+T1·(1-r+r²)) = R4",
   fabs(R_from_S2 - R4) < mpf('1e-50'))

# A60: R∞ gap < R3 gap (infinite sum is NOT closer to TARGET than R3)
ck("A60: R∞ gap is comparable to R3 gap",
   gap_inf < gap_R3 * 2)

# ─────────────────────────────────────────────────────────────────────────────
# Section 6: Additional algebraic identities
# ─────────────────────────────────────────────────────────────────────────────

# A61: r·μ² = h·E_μ  (defining relation, exact)
ck("A61: r·μ² = h·E_μ",
   fabs(r * ALPHA_INV**2 - h*E_mu) < mpf('1e-58'))

# A62: r·μ² = h·e²  (since E_μ = e²)
ck("A62: r·μ² = h·e²",
   fabs(r * ALPHA_INV**2 - h*E_e**2) < mpf('1e-58'))

# A63: (1+r)⁻¹ = μ²/(μ²+h·E_μ)  (the resummation denominator)
ck("A63: (1+r)⁻¹ = μ²/(μ²+hE_μ)",
   fabs(1/(1+r) - ALPHA_INV**2/(ALPHA_INV**2+h*E_mu)) < mpf('1e-55'))

# A64: T1·(1+r)⁻¹ = hμ/(e³(μ²+hE_μ))
lhs = T1/(1+r)
rhs = h*ALPHA_INV/(E_e**3*(ALPHA_INV**2+h*E_mu))
ck("A64: T1/(1+r) = hμ/(e³(μ²+hE_μ))",
   fabs(lhs - rhs) < mpf('1e-55'))

# A65: r = (μ²+hEμ - μ²)/μ² = Δμ/μ²  (r measures the fractional monad deformation)
delta_mu = h*E_mu
ck("A65: r = h·E_μ/μ²  (fractional deformation of the monad squared)",
   fabs(r - delta_mu/ALPHA_INV**2) < mpf('1e-58'))

# ─────────────────────────────────────────────────────────────────────────────
# Final summary
# ─────────────────────────────────────────────────────────────────────────────
print()
print("=== P214 Key Results ===")
print()
print(f"  h∨(A₂) = {h}   (dual Coxeter number of A₂)")
print(f"  E_e = π = {nstr(E_e, 10)}")
print(f"  E_μ = π² = {nstr(E_mu, 10)}")
print(f"  μ   = {nstr(ALPHA_INV, 15)}")
print(f"  α   = {nstr(alpha, 15)}")
print()
print(f"  r = h·E_μ·α² = {nstr(r, 25)}")
print(f"    = 3π²/μ²    = {nstr(3*pi**2/ALPHA_INV**2, 25)}")
print()
print("  Proposition P214.1 (PROVEN):")
print(f"    T₂/T₁ = -r: diff = {nstr(fabs(T2/T1 - (-r)), 5)}")
print(f"    T₃/T₂ = -r: diff = {nstr(fabs(T3/T2 - (-r)), 5)}")
print(f"    T₄/T₃ = -r: diff = {nstr(fabs(T4/T3 - (-r)), 5)}")
print()
print("  Conjecture P214.2 (numerically verified):")
print(f"    h∨(A₂)·E_μ·α² = r: diff = {nstr(fabs(h*E_mu*alpha**2 - r), 5)}")
print(f"    (3 unoriented 2-step paths on A₂ triangle, each E_e²·α²)")
print()
print(f"  R∞ = {nstr(R_inf, 22)}")
print(f"  Gap = {nstr(gap_inf, 10)} ppm")
print(f"  Deformed monad denominator: μ²+h·E_μ = μ²(1+r) = {nstr(ALPHA_INV**2*(1+r), 15)}")
print()
print("  Status: OP-G2A2-ratio CONJECTURED")
print("  Algebraic derivation complete (Prop. P214.1).")
print("  Geometric origin of ±2 exponent shifts: open.")

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