"""
verify_P162.py — Verification script for Addendum 162 (correction to P87)

Checks:
  1. Each x_i^2 * w_i term with correct exponents
  2. M2_correct ≈ 13929 (within 1%)
  3. M2_paper  ≈ 354    (the erroneous value from P87 Prop 2.1)
  4. ratio M2_correct / M2_paper ≈ 39

Copyright: Léon Fernando Vlegels. License: MIT. May 2026.
"""

import sys

import mpmath

mpmath.mp.dps = 50

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

pi    = mpmath.pi
mu0   = 4*pi**3 + pi**2 + pi          # = α⁻¹ ≈ 137.036  (verified against P87 §1)

# ── Eigenvalues and weights ────────────────────────────────────────────────
x1, x2, x3 = pi, pi**2, 4*pi**3       # sector eigenvalues
w1, w2, w3 = x1/mu0, x2/mu0, x3/mu0  # spectral weights  w_i = x_i / μ₀

print("=== P162 Verification ===\n")
print(f"μ₀  = {float(mu0):.6f}  (expect ≈ 137.036)")
print(f"x₁  = π    = {float(x1):.6f}")
print(f"x₂  = π²   = {float(x2):.6f}")
print(f"x₃  = 4π³  = {float(x3):.6f}")
print(f"w₁  = {float(w1):.6f}")
print(f"w₂  = {float(w2):.6f}")
print(f"w₃  = {float(w3):.6f}")
print(f"Σwᵢ = {float(w1+w2+w3):.10f}  (must be 1)")
print()

# ── Correct term-by-term computation ─────────────────────────────────────
term1_correct = x1**2 * w1   # = π² · (π/μ₀)  = π³/μ₀
term2_correct = x2**2 * w2   # = π⁴ · (π²/μ₀) = π⁶/μ₀
term3_correct = x3**2 * w3   # = 16π⁶ · (4π³/μ₀) = 64π⁹/μ₀

M2_correct = term1_correct + term2_correct + term3_correct

print("── Correct exponents ─────────────────────────────────────────────────")
print(f"x₁²·w₁ = π³/μ₀        = {float(term1_correct):.6f}  (expect π³/μ₀ ≈ {float(pi**3/mu0):.4f})")
print(f"x₂²·w₂ = π⁶/μ₀        = {float(term2_correct):.6f}  (expect π⁶/μ₀ ≈ {float(pi**6/mu0):.4f})")
print(f"x₃²·w₃ = 64π⁹/μ₀      = {float(term3_correct):.4f}  (expect 64π⁹/μ₀ ≈ {float(64*pi**9/mu0):.2f})")
print(f"M₂_correct = {float(M2_correct):.4f}")
print()

# ── P87 Proposition 2.1 — erroneous formula ───────────────────────────────
# Wrong: x₂²w₂ treated as π⁵/μ₀ (exponent 5 instead of 6)
#        x₃²w₃ treated as 16π⁷/μ₀ (coefficient 16, exponent 7 instead of 64π⁹)
term1_paper = pi**3 / mu0
term2_paper = pi**5 / mu0   # ← wrong: should be π⁶/μ₀
term3_paper = 16*pi**7 / mu0  # ← wrong: should be 64π⁹/μ₀

M2_paper = term1_paper + term2_paper + term3_paper

print("── P87 Prop 2.1 (erroneous) ──────────────────────────────────────────")
print(f"x₁²·w₁ = π³/μ₀        = {float(term1_paper):.6f}")
print(f"x₂²·w₂ = π⁵/μ₀ (WRONG)= {float(term2_paper):.6f}  (correct: {float(term2_correct):.4f})")
print(f"x₃²·w₃ = 16π⁷/μ₀ (WRONG)= {float(term3_paper):.4f}  (correct: {float(term3_correct):.2f})")
print(f"M₂_paper = {float(M2_paper):.4f}")
print()

# ── Assertions ────────────────────────────────────────────────────────────
ratio = M2_correct / M2_paper

print("S1  Assertions")

# 1. Correct terms match closed-form powers of π
check(1, "term2 = π⁶/μ₀  and  term3 = 64π⁹/μ₀  (exact)",
      abs(float(term2_correct - pi**6/mu0)) < 1e-40
      and abs(float(term3_correct - 64*pi**9/mu0)) < 1e-40)

# 2. M2_correct ≈ 13929 within 1%
expected_correct = 13929
check(2, f"M₂_correct = {float(M2_correct):.1f} ≈ {expected_correct} (within 1%)",
      abs(float(M2_correct) - expected_correct) / expected_correct < 0.01)

# 3. M2_paper ≈ 354 within 1%  (the wrong value)
expected_paper = 354
check(3, f"M₂_paper  = {float(M2_paper):.1f} ≈ {expected_paper} (within 1%)",
      abs(float(M2_paper) - expected_paper) / expected_paper < 0.01)

# 4. Ratio ≈ 39
expected_ratio = 39
check(4, f"M₂_correct/M₂_paper = {float(ratio):.2f} ≈ {expected_ratio} (factor-of-39 error confirmed)",
      abs(float(ratio) - expected_ratio) / expected_ratio < 0.05)

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