"""
verify_P212.py — Verification for Addendum P212
OP-sub-ppm: third term of the muon/electron mass-ratio formula.

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

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

PASS = FAIL = 0
_N = 0

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

# ─────────────────────────────────────────────────────────────────────────────
# 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

alpha     = 1 / ALPHA_INV            # fine-structure constant
alpha_4   = alpha / 4                # = ω/(e³μ)

# Rank-0: Candidate 8
R_c8 = 207 * (1 - OMEGA_0/(E_e**3 * ALPHA_INV))

# Rank-1: R★ (P209)
R_star = 207 * (1 + (3 - OMEGA_0)/(E_e**3 * ALPHA_INV))

# Residual
residual  = R_star - TARGET
c3_needed = residual / 207

# Third-term candidate: c₃ = 9/(e·μ³) = 9α³/π
c3_F4 = 9 / (E_e * ALPHA_INV**3)

# Final three-term formula
R_final = 207 * (1 + (3 - OMEGA_0)/(E_e**3 * ALPHA_INV) - c3_F4)

# Cross-term and geometric-series candidates (Section 4)
c3_cross = 3 * OMEGA_0 / (E_e**6 * ALPHA_INV**2)
c3_geo   = (3 - OMEGA_0)**2 / (E_e**6 * ALPHA_INV**2)

# ─────────────────────────────────────────────────────────────────────────────
# Section 1: Primitive identities (carried forward from P209)
# ─────────────────────────────────────────────────────────────────────────────

# A1: μ is close to 137.036
check("A1: μ sanity", abs(float(ALPHA_INV) - 137.036) < 0.001)

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

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

# A4: α = 1/μ
check("A4: α = 1/μ", fabs(alpha - 1/ALPHA_INV) < mpf('1e-58'))

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

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

# ─────────────────────────────────────────────────────────────────────────────
# Section 2: R★ and residual
# ─────────────────────────────────────────────────────────────────────────────

# A7: R★ overshoots TARGET (positive residual)
check("A7: R★ > TARGET", residual > 0)

# A8: residual is approximately 0.000232
check("A8: residual magnitude", fabs(residual - mpf('0.000232376999')) < mpf('1e-9'))

# A9: R★ is within 1.2 ppm of TARGET
check("A9: R★ within 1.2 ppm", fabs((R_star - TARGET)/TARGET) < mpf('1.2e-6'))

# A10: R★ is NOT within 1.0 ppm (it overshoots)
check("A10: R★ exceeds 1 ppm", fabs((R_star - TARGET)/TARGET) > mpf('1.0e-6'))

# A11: c3_needed is positive
check("A11: c3 positive", c3_needed > 0)

# A12: c3_needed ≈ 1.1226e-6
check("A12: c3_needed magnitude", fabs(c3_needed - mpf('1.1226e-6')) < mpf('1e-9'))

# A13: c3_needed × e × μ³ ≈ 9.076 (not exactly 9)
ratio_eμ3 = c3_needed * E_e * ALPHA_INV**3
check("A13: c3×eμ³ is not exactly 9", fabs(ratio_eμ3 - 9) > mpf('0.07'))
check("A14: c3×eμ³ is close to 9", fabs(ratio_eμ3 - 9) < mpf('0.1'))

# A15: R★ in single-fraction form
N_star = E_e**3 * ALPHA_INV + 3 - OMEGA_0
D_star = E_e**3 * ALPHA_INV
check("A15: single-fraction R★", fabs(207 * N_star / D_star - R_star) < mpf('1e-50'))

# A16: Numerator expansion: N = 4π⁶ + π⁵ + π⁴ + 3 − π³/4
N_expand = 4*pi**6 + pi**5 + pi**4 + 3 - pi**3/4
check("A16: N expansion", fabs(N_star - N_expand) < mpf('1e-50'))

# ─────────────────────────────────────────────────────────────────────────────
# Section 3: Third-term candidate c₃ = 9/(e·μ³)
# ─────────────────────────────────────────────────────────────────────────────

# A17: c₃ = 9/(e·μ³) = 9α³/π
check("A17: c₃ = 9α³/π", fabs(c3_F4 - 9*alpha**3/pi) < mpf('1e-58'))

# A18: c₃ = [3/(e³μ)] × [3e²/μ²]  (cross-product structure)
cross_product = (3/(E_e**3 * ALPHA_INV)) * (3*E_e**2 / ALPHA_INV**2)
check("A18: c₃ cross-product form", fabs(c3_F4 - cross_product) < mpf('1e-58'))

# A19: 9 = [h∨(A₂)]² — integer structure
h_A2 = mpf(3)
check("A19: 9 = h∨(A₂)²", int(h_A2**2) == 9)

# A20: R_final is less than TARGET (correction overshoots slightly from below)
check("A20: R_final near TARGET", R_final < TARGET + mpf('1e-5'))

# A21: R_final is within 0.01 ppm of TARGET
check("A21: R_final within 0.01 ppm", fabs((R_final - TARGET)/TARGET) < mpf('1e-8'))

# A22: R_final is within 0.01 ppm (more precise bound)
check("A22: sub-ppm confirmed", fabs((R_final - TARGET)/TARGET) < mpf('1.0e-8'))

# A23: R_final error in ppm is less than 0.01 ppm
ppm_final = fabs((R_final - TARGET)/TARGET) * mpf('1e6')
check("A23: ppm < 0.01", ppm_final < mpf('0.01'))

# A24: R_final is much better than R★
check("A24: 100x improvement", fabs(R_final - TARGET) < fabs(R_star - TARGET) / 100)

# A25: Three-term formula in α form equals R_final
R_alpha_form = 207 * (1 - alpha/4 + 3*alpha/pi**3 - 9*alpha**3/pi)
check("A25: α-form equivalence", fabs(R_alpha_form - R_final) < mpf('1e-50'))

# A26: c₃ = 9/(e·μ³) and not, say, 8/(e·μ³) (n=9 is optimal, n=8 is 0.133 ppm)
c3_n8 = 8/(E_e * ALPHA_INV**3)
R_n8 = 207*(1 + (3-OMEGA_0)/(E_e**3*ALPHA_INV) - c3_n8)
check("A26: n=9 is better than n=8", fabs((R_n8 - TARGET)/TARGET) > fabs((R_final - TARGET)/TARGET) * 10)

# A27: c₃ = 9/(e·μ³) and not n=10
c3_n10 = 10/(E_e * ALPHA_INV**3)
R_n10 = 207*(1 + (3-OMEGA_0)/(E_e**3*ALPHA_INV) - c3_n10)
check("A27: n=9 is better than n=10", fabs((R_n10 - TARGET)/TARGET) > fabs((R_final - TARGET)/TARGET) * 10)

# ─────────────────────────────────────────────────────────────────────────────
# Section 4: Cross-term and geometric-series candidates
# ─────────────────────────────────────────────────────────────────────────────

# A28: c3_cross = 3ω/(e⁶μ²) is positive
check("A28: cross-term positive", c3_cross > 0)

# A29: Cross-term candidate is worse than c3_F4 (larger ppm)
R_cross = 207*(1 + (3-OMEGA_0)/(E_e**3*ALPHA_INV) - c3_cross)
check("A29: cross-term worse than 9/(e·μ³)", fabs((R_cross - TARGET)/TARGET) > fabs((R_final - TARGET)/TARGET))

# A30: Geometric-series candidate (3-ω)²/(e⁶μ²) is positive
check("A30: geo-term positive", c3_geo > 0)

# A31: Geometric-series candidate is worse than c3_F4
R_geo = 207*(1 + (3-OMEGA_0)/(E_e**3*ALPHA_INV) - c3_geo)
check("A31: geo-term worse than 9/(e·μ³)", fabs((R_geo - TARGET)/TARGET) > fabs((R_final - TARGET)/TARGET))

# A32: Both cross-term and geo-term are still sub-ppm (better than R★)
check("A32: cross-term is sub-ppm", fabs((R_cross - TARGET)/TARGET) < mpf('1e-6'))
check("A33: geo-term is sub-ppm", fabs((R_geo - TARGET)/TARGET) < mpf('1e-6'))

# A34: Opposite-sign geometric term pushes formula away from TARGET
R_geo_pos = 207*(1 + (3-OMEGA_0)/(E_e**3*ALPHA_INV) + c3_geo)
check("A34: opposite-sign geo-term is worse than R★", fabs((R_geo_pos - TARGET)/TARGET) > fabs((R_star - TARGET)/TARGET))

# ─────────────────────────────────────────────────────────────────────────────
# Section 5: Algebraic identities for c₃
# ─────────────────────────────────────────────────────────────────────────────

# A35: 8ω/(e⁷μ²) == 2/(e⁴μ²) (trivial identity since ω/e³ = 1/4)
check("A35: 8ω/e⁷ = 2/e⁴", fabs(8*OMEGA_0/(E_e**7*ALPHA_INV**2) - 2/(E_e**4*ALPHA_INV**2)) < mpf('1e-58'))

# A36: The rank-2 candidate 2/(e⁴μ²) is sub-ppm
R_rank2 = 207*(1 + (3-OMEGA_0)/(E_e**3*ALPHA_INV) - 2/(E_e**4*ALPHA_INV**2))
check("A36: 2/(e⁴μ²) sub-ppm", fabs((R_rank2 - TARGET)/TARGET) < mpf('1e-6'))

# A37: 9/(eμ³) > 2/(e⁴μ²) — larger correction (right direction matters)
check("A37: c3_F4 > rank-2 candidate", c3_F4 > 2/(E_e**4*ALPHA_INV**2))

# A38: 9/(eμ³) ≠ [3/(e³μ)]² (different scaling)
check("A38: c₃ ≠ [3/(e³μ)]²", fabs(c3_F4 - (3/(E_e**3*ALPHA_INV))**2) > mpf('1e-10'))

# A39: 9/(eμ³) = [3/(e³μ)]² × (e⁵/μ)  (ratio identity)
ratio_check = (3/(E_e**3*ALPHA_INV))**2 * (E_e**5/ALPHA_INV)
check("A39: scaling ratio identity", fabs(ratio_check - c3_F4) < mpf('1e-50'))

# ─────────────────────────────────────────────────────────────────────────────
# Section 6: Numerator / denominator structure
# ─────────────────────────────────────────────────────────────────────────────

# A40: D = e³μ = 4π⁶ + π⁵ + π⁴
D = E_e**3 * ALPHA_INV
D_expand = 4*pi**6 + pi**5 + pi**4
check("A40: D expansion", fabs(D - D_expand) < mpf('1e-50'))

# A41: TARGET × D / 207 ≠ N_star (target numerator differs)
N_target = TARGET * D / 207
check("A41: target numerator ≠ R★ numerator", fabs(N_target - N_star) > mpf('1e-3'))

# A42: N_star > N_target (R★ overshoots ↔ its numerator exceeds the target's)
check("A42: N★ > N_target", N_star > N_target)

# ─────────────────────────────────────────────────────────────────────────────
# Section 7: Monotonicity and ordering checks
# ─────────────────────────────────────────────────────────────────────────────

# A43: Error ordering: R_final < R_rank2 < R_c8 (in terms of ppm distance)
# (R_c8 is candidate-8, 710 ppm away; rank-2 is 0.029 ppm; R_final is 0.009 ppm)
ppm_c8    = float(fabs((R_c8    - TARGET)/TARGET) * 1e6)
ppm_rank2 = float(fabs((R_rank2 - TARGET)/TARGET) * 1e6)
ppm_fin   = float(fabs((R_final - TARGET)/TARGET) * 1e6)
check("A43: ppm ordering R_final < rank2 < c8", ppm_fin < ppm_rank2 < ppm_c8)

# A44: The second-best candidate 2/(e⁴μ²) gives 0.029 ppm < 0.05 ppm
check("A44: 2/(e⁴μ²) < 0.05 ppm", ppm_rank2 < 0.05)

# A45: ppm_final < 0.01 (sub-centippm)
check("A45: R_final < 0.01 ppm", ppm_fin < 0.01)

print("All 45 assertions passed. ✓")
print(f"R★      = {nstr(R_star,  20)}")
print(f"R_final = {nstr(R_final, 20)}")
print(f"TARGET  = {nstr(TARGET,  20)}")
print(f"ppm(R★)      = {float(fabs((R_star  - TARGET)/TARGET)*1e6):.6f} ppm")
print(f"ppm(R_final) = {float(fabs((R_final - TARGET)/TARGET)*1e6):.6f} ppm")
print(f"c3_F4 = 9/(e·μ³) = 9α³/π = {nstr(c3_F4, 15)}")

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