"""
verify_P213.py — Verification for Addendum P213
G₂–A₂ Unified Mass-Ratio Theorem: geometric series, closed form, fourth-term probe.

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

import sys

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

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}")

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

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

# A1: G₂ dimension check: dim(G₂) = 14
check(1, "A1: 14·15 − h∨(A₂) = 207", 14*15 - int(h) == 207)

# A2: 207 = 14·15 − 3
check(2, "A2: 207 = 14·15 − 3", 14*15 - 3 == 207)

# A3: h∨(A₂) = 3
check(3, "A3: h∨(A₂) = 3", int(h) == 3)

# A4: μ ≈ 137.036
check(4, "A4: μ ≈ 137.036", fabs(ALPHA_INV - mpf('137.036')) < mpf('0.001'))

# A5: ω = π³/4
check(5, "A5: ω = π³/4", fabs(OMEGA_0 - E_e**3/4) < mpf('1e-58'))

# A6: ω/e³ = 1/4
check(6, "A6: ω/e³ = 1/4", fabs(OMEGA_0/E_e**3 - mpf('1')/4) < mpf('1e-58'))

# A7: μ = 4e³ + e² + e
check(7, "A7: μ = 4e³+e²+e", fabs(ALPHA_INV - (4*E_e**3 + E_e**2 + E_e)) < mpf('1e-58'))

# A8: α·μ = 1
check(8, "A8: α·μ = 1", fabs(alpha * ALPHA_INV - 1) < mpf('1e-58'))

# ─────────────────────────────────────────────────────────────────────────────
# Section 1: Series terms T_k = (-1)^{k+1} h^k e^{2k-5} / μ^{2k-1}
# ─────────────────────────────────────────────────────────────────────────────
T0 = -OMEGA_0 / (E_e**3 * ALPHA_INV)          # −ω/e³μ  (candidate-8 base)
T1 = +h      / (E_e**3 * ALPHA_INV)           # +h/e³μ
T2 = -h**2   / (E_e    * ALPHA_INV**3)        # −h²/eμ³
T3_pred = +h**3 * E_e / ALPHA_INV**5          # +h³e/μ⁵  (geometric prediction)

# A9: T0 = −1/(4μ)
check(9, "A9: T0 = −1/4μ", fabs(T0 - (-1/(4*ALPHA_INV))) < mpf('1e-58'))

# A10: T0 = −ω/e³μ = −α/4
check(10, "A10: T0 = −α/4", fabs(T0 - (-alpha/4)) < mpf('1e-58'))

# A11: T1 = h/(e³μ) = 3α/π³ = 3α/e³
check(11, "A11: T1 = 3α/e³", fabs(T1 - 3*alpha/E_e**3) < mpf('1e-58'))

# A12: T2 = −h²/(eμ³) = −9α³/π = −9α³/e
check(12, "A12: T2 = −9α³/e", fabs(T2 - (-9*alpha**3/E_e)) < mpf('1e-55'))

# A13: T3_pred = h³·e/μ⁵ = 27π·α⁵/π = 27α⁵/e⁴... let's verify directly
# T3_pred = 27·e/μ⁵ = 27·α⁵·e/α⁵... = 27 * E_e / ALPHA_INV**5
check(13, "A13: T3_pred = 27e/μ⁵", fabs(T3_pred - 27*E_e/ALPHA_INV**5) < mpf('1e-58'))

# A14: T1 is positive
check(14, "A14: T1 > 0", T1 > 0)

# A15: T2 is negative
check(15, "A15: T2 < 0", T2 < 0)

# A16: T3_pred is positive
check(16, "A16: T3_pred > 0", T3_pred > 0)

# A17: |T1| > |T2| > |T3_pred|  (series is monotone decreasing in magnitude)
check(17, "A17: |T1| > |T2| > |T3_pred|", fabs(T1) > fabs(T2) > fabs(T3_pred))

# ─────────────────────────────────────────────────────────────────────────────
# Section 2: Geometric ratio r = h·e²/μ² = 3π²/μ²
# ─────────────────────────────────────────────────────────────────────────────
r = h * E_e**2 / ALPHA_INV**2    # = 3π²/ALPHA_INV²

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

# A19: r is small — roughly 0.00158
check(19, "A19: r ≈ 0.001577", fabs(r - mpf('0.001577')) < mpf('0.00001'))

# A20: |T2/T1| = r exactly (geometric ratio at step k=2/1)
check(20, "A20: |T2/T1| = r", fabs(fabs(T2/T1) - r) < mpf('1e-55'))

# A21: |T3_pred/T2| = r exactly (geometric ratio at step k=3/2)
check(21, "A21: |T3_pred/T2| = r", fabs(fabs(T3_pred/T2) - r) < mpf('1e-55'))

# A22: the ratio is the SAME at both steps — geometric series confirmed
check(22, "A22: ratio identical at both steps", fabs(fabs(T2/T1) - fabs(T3_pred/T2)) < mpf('1e-60'))

# A23: T2/T1 is negative (sign alternates)
check(23, "A23: T2/T1 < 0 (sign alternates)", T2/T1 < 0)

# A24: T3_pred/T2 is negative (sign alternates again)
check(24, "A24: T3_pred/T2 < 0 (sign alternates)", T3_pred/T2 < 0)

# A25: common ratio T2/T1 = T3_pred/T2 = −r
check(25, "A25: T2/T1 = −r", fabs(T2/T1 - (-r)) < mpf('1e-55'))
check(26, "A26: T3_pred/T2 = −r", fabs(T3_pred/T2 - (-r)) < mpf('1e-55'))

# ─────────────────────────────────────────────────────────────────────────────
# Section 3: R-value computations
# ─────────────────────────────────────────────────────────────────────────────
R_c8    = 207*(1 + T0)
R_star  = 207*(1 + T0 + T1)
R3      = 207*(1 + T0 + T1 + T2)
R4_pred = 207*(1 + T0 + T1 + T2 + T3_pred)

gap_c8    = fabs(R_c8    - TARGET)/TARGET * 1000000   # in ppm
gap_star  = fabs(R_star  - TARGET)/TARGET * 1000000
gap_R3    = fabs(R3      - TARGET)/TARGET * 1000000
gap_R4    = fabs(R4_pred - TARGET)/TARGET * 1000000

# A27: 14·15 − 3 = 207 in the prefactor
check(27, "A27: prefactor = 14·15−h∨(A₂)", fabs(R_c8 - (14*15 - 3)*(1 + T0)) < mpf('1e-50'))

# A28: R_c8 gap is ~706 ppm (≈ 0.071%)
check(28, "A28: R_c8 gap ≈ 706 ppm", gap_c8 > mpf('700') and gap_c8 < mpf('710'))

# A29: R_c8 < TARGET (candidate-8 undershoots)
check(29, "A29: R_c8 < TARGET", R_c8 < TARGET)

# A30: R_star > TARGET (R★ overshoots)
check(30, "A30: R_star > TARGET", R_star > TARGET)

# A31: R_star gap is in [1.1, 1.2] ppm
check(31, "A31: R_star gap ∈ [1.1, 1.2] ppm", gap_star > mpf('1.1') and gap_star < mpf('1.2'))

# A32: R3 > TARGET (three-term formula still slightly overshoots)
check(32, "A32: R3 > TARGET", R3 > TARGET)

# A33: R3 gap < 0.01 ppm
check(33, "A33: R3 gap < 0.01 ppm", gap_R3 < mpf('0.01'))

# A34: R3 gap is approximately 0.00937 ppm
check(34, "A34: R3 gap ≈ 0.009373 ppm", fabs(gap_R3 - mpf('0.009373')) < mpf('0.0001'))

# A35: R3 is 100× closer to TARGET than R_star
check(35, "A35: R3 100× better than R_star", fabs(R3 - TARGET) < fabs(R_star - TARGET) / 100)

# A36: R3 is 75000× closer to TARGET than R_c8
check(36, "A36: R3 >> R_c8", fabs(R3 - TARGET) < fabs(R_c8 - TARGET) / 70000)

# A37: The predicted T3 makes R4 WORSE than R3 (series overshoots)
check(37, "A37: R4_pred gap > R3 gap (T3_pred overshoots)", gap_R4 > gap_R3)

# A38: R4_pred > R3 > TARGET
check(38, "A38: R4_pred > R3 > TARGET", R4_pred > R3 > TARGET)

# A39: R4_pred ≈ R_closed (differs by series convergence order r²)
# R_closed = 207*(1 + T0 + T1/(1+r)) — check via direct formula
R_closed = 207*(1 + T0 + T1/(1+r))
check(39, "A39: R4_pred ≈ R_closed", fabs(R4_pred - R_closed) < fabs(T3_pred) * 207 * r)

# ─────────────────────────────────────────────────────────────────────────────
# Section 4: Closed-form geometric sum
# ─────────────────────────────────────────────────────────────────────────────
S_closed = T1 / (1 + r)
R_closed = 207*(1 + T0 + S_closed)

gap_closed = fabs(R_closed - TARGET)/TARGET * 1000000

# A40: S_closed = T1·μ²/(μ²+h·e²) (explicit form)
S_explicit = T1 * ALPHA_INV**2 / (ALPHA_INV**2 + h*E_e**2)
check(40, "A40: S_closed explicit form", fabs(S_closed - S_explicit) < mpf('1e-58'))

# A41: S_closed = h·μ/(e³(μ²+he²))
S_form2 = h * ALPHA_INV / (E_e**3 * (ALPHA_INV**2 + h*E_e**2))
check(41, "A41: S_closed = hμ/(e³(μ²+he²))", fabs(S_closed - S_form2) < mpf('1e-55'))

# A42: Fully explicit closed form
S_form3 = 3*ALPHA_INV / (E_e**3 * (ALPHA_INV**2 + 3*E_e**2))
check(42, "A42: S_closed = 3μ/(e³(μ²+3e²))", fabs(S_closed - S_form3) < mpf('1e-55'))

# A43: R_closed = 207·(1 − ω/e³μ + hμ/(e³(μ²+he²)))
R_closed_check = 207*(1 - OMEGA_0/(E_e**3*ALPHA_INV) + h*ALPHA_INV/(E_e**3*(ALPHA_INV**2+h*E_e**2)))
check(43, "A43: R_closed explicit form", fabs(R_closed - R_closed_check) < mpf('1e-50'))

# A44: R_closed > TARGET (the infinite sum also overshoots)
check(44, "A44: R_closed > TARGET", R_closed > TARGET)

# A45: R_closed gap ≈ 0.01113 ppm
check(45, "A45: R_closed gap ≈ 0.01113 ppm", fabs(gap_closed - mpf('0.01113')) < mpf('0.0001'))

# A46: R_closed > R3 (infinite sum slightly worse than N=2 truncation)
check(46, "A46: R_closed > R3", R_closed > R3)

# A47: R_closed − R3 = 207·T1·r²/(1+r)  (geometric series partial-sum error formula)
delta_cr = R_closed - R3
expected_delta = 207 * T1 * r**2 / (1+r)
check(47, "A47: R_closed − R3 = 207·T1·r²/(1+r)", fabs(delta_cr - expected_delta) < mpf('1e-50'))

# A48: |r| < 0.002 (series convergence — ratio well below 1)
check(48, "A48: r < 0.002", r < mpf('0.002'))

# A49: S_closed/T1 < 1.002 (geometric series sums to near T1)
check(49, "A49: S_closed ≈ T1 to within 2r", fabs(S_closed/T1 - 1) < r * 2)

# ─────────────────────────────────────────────────────────────────────────────
# Section 5: Fourth-term grid scan — best T3 candidates
# ─────────────────────────────────────────────────────────────────────────────
# Best from grid: T3_best = −h³/e³/μ⁴ (sign=−1, a=−3, b=4)
T3_neg_e3_mu4 = -h**3 / (E_e**3 * ALPHA_INV**4)   # −27/(e³μ⁴)
T3_neg_e1_mu5 = -h**3 * E_e / ALPHA_INV**5         # −27e/μ⁵ = −T3_pred
T3_neg_mu5    = -h**3 / ALPHA_INV**5               # −27/μ⁵

R4_grid1 = 207*(1 + T0 + T1 + T2 + T3_neg_e3_mu4)
R4_grid2 = 207*(1 + T0 + T1 + T2 + T3_neg_e1_mu5)
R4_grid3 = 207*(1 + T0 + T1 + T2 + T3_neg_mu5)

gap_grid1 = fabs(R4_grid1 - TARGET)/TARGET * 1e6
gap_grid2 = fabs(R4_grid2 - TARGET)/TARGET * 1e6
gap_grid3 = fabs(R4_grid3 - TARGET)/TARGET * 1e6

# A50: Best grid candidate improves on R3
check(50, "A50: −27/(e³μ⁴) improves on R3", gap_grid1 < gap_R3)

# A51: Best candidate gap < 0.007 ppm
check(51, "A51: best T3 gap < 0.007 ppm", gap_grid1 < mpf('0.0070'))

# A52: The negative T3_pred is second-best
check(52, "A52: −T3_pred also improves on R3", gap_grid2 < gap_R3)

# A53: T3_neg_e3_mu4 outperforms T3_neg_e1_mu5
check(53, "A53: −27/(e³μ⁴) < −27e/μ⁵", gap_grid1 < gap_grid2)

# A54: The positive predicted T3 is WORSE than R3
check(54, "A54: +T3_pred worsens R3", gap_R4 > gap_R3)

# A55: The sign-alternating series with T3_pred goes in the wrong direction
# (because R3 already overshoots TARGET; T3_pred is positive and pushes further)
check(55, "A55: sign-alternating T3 over-corrects beyond R3", R4_pred > R3 > TARGET)

# ─────────────────────────────────────────────────────────────────────────────
# Section 6: Residual analysis
# ─────────────────────────────────────────────────────────────────────────────
Delta3 = R3 - TARGET

# A56: Δ₃ is positive and small
check(56, "A56: Δ₃ > 0", Delta3 > 0)
check(57, "A57: Δ₃ < 2×10⁻⁶", Delta3 < mpf('2e-6'))

# A58: Δ₃ ≈ 1.938×10⁻⁶
check(58, "A58: Δ₃ ≈ 1.938×10⁻⁶", fabs(Delta3 - mpf('1.938e-6')) < mpf('1e-9'))

# A59: Δ₃ / |T2| ≈ 1.741  (residual is ~1.74× the third term)
ratio_delta_T2 = Delta3 / fabs(T2)
check(59, "A59: Δ₃/|T2| ≈ 1.741", fabs(ratio_delta_T2 - mpf('1.741')) < mpf('0.001'))

# A60: Δ₃ is NOT equal to T3_pred × 207  (geometric prediction is off)
check(60, "A60: Δ₃ ≠ 207·T3_pred", fabs(Delta3 - 207*T3_pred) > mpf('1e-7'))

# ─────────────────────────────────────────────────────────────────────────────
# Section 7: Algebraic identities for the ratio r
# ─────────────────────────────────────────────────────────────────────────────

# A61: r = 3π²/μ² = 3α²π² = h·α²·e²
check(61, "A61: r = h·α²·e²", fabs(r - h*alpha**2*E_e**2) < mpf('1e-55'))

# A62: r = T2/T1 × (−1) — sign-stripped ratio
check(62, "A62: r = |T2/T1|", fabs(r - fabs(T2/T1)) < mpf('1e-55'))

# A63: 1 + r ≈ 1 to within 0.16%
check(63, "A63: r is small", fabs(1 + r - 1) < mpf('0.002'))

# A64: The closed-form denominator μ²+he² = μ²(1+r)
denom = ALPHA_INV**2 + h*E_e**2
check(64, "A64: μ²+he² = μ²(1+r)", fabs(denom - ALPHA_INV**2*(1+r)) < mpf('1e-50'))

# A65: Partial sum at N=2 of geometric series T1·Σ(−r)^k for k=0,1,2:
S2 = T1*(1 - r + r**2)     # geometric partial sum order 2
R_partsum2 = 207*(1 + T0 + S2)
check(65, "A65: partial sum N=2 = R4_pred", fabs(R_partsum2 - R4_pred) < mpf('1e-50'))

# ─────────────────────────────────────────────────────────────────────────────
# Final summary
# ─────────────────────────────────────────────────────────────────────────────
if FAIL == 0:
    print("All 65 assertions passed. ✓")
print()
print(f"  ALPHA_INV (μ) = {nstr(ALPHA_INV, 20)}")
print(f"  h = h∨(A₂)   = {h}")
print(f"  r = h·e²/μ²  = {nstr(r, 20)}")
print()
print(f"  R_c8    = {nstr(R_c8,    22)}   gap = {nstr(gap_c8,    8)} ppm")
print(f"  R_star  = {nstr(R_star,  22)}   gap = {nstr(gap_star,  8)} ppm")
print(f"  R3      = {nstr(R3,      22)}   gap = {nstr(gap_R3,    8)} ppm")
print(f"  R4_pred = {nstr(R4_pred, 22)}   gap = {nstr(gap_R4,    8)} ppm  [+T3_pred, WORSE]")
print(f"  R_closed= {nstr(R_closed,22)}   gap = {nstr(gap_closed,8)} ppm  [infinite sum]")
print(f"  R4_best = {nstr(R4_grid1,22)}   gap = {nstr(gap_grid1, 8)} ppm  [-27/e³μ⁴]")
print()
print(f"  Geometric ratio r = {nstr(r, 15)}")
print(f"  |T2/T1| − r = {nstr(fabs(fabs(T2/T1)-r), 5)}")
print(f"  |T3/T2| − r = {nstr(fabs(fabs(T3_pred/T2)-r), 5)}")
print()
print(f"  R3 residual Δ₃ = {nstr(R3 - TARGET, 15)}")

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