"""
verify_P087.py — P87: The Second Spectral Moment and Tr(X_sector²)

Paper claims:  M₂ = (π³+π⁵+16π⁷)/μ₀ ≈ 353.86

GENUINE PAPER ERROR DETECTED:
The discrete second moment is M₂ = Σ x_i²·w_i = Σ x_i³/μ₀
with x₁=π, x₂=π², x₃=4π³ and w_i=x_i/μ₀.

Correct computation:
  x₁²·w₁ = π²·(π/μ₀) = π³/μ₀           ✓ (paper has this right)
  x₂²·w₂ = (π²)²·(π²/μ₀) = π⁶/μ₀      but paper writes π⁵/μ₀  (WRONG: 4+2=6, not 5)
  x₃²·w₃ = (4π³)²·(4π³/μ₀) = 64π⁹/μ₀  but paper writes 16π⁷/μ₀ (WRONG)

Correct: M₂ = (π³+π⁶+64π⁹)/μ₀ ≈ 13929
Paper:   M₂ = (π³+π⁵+16π⁷)/μ₀ ≈ 355  (arithmetic error in exponent addition)

Additional paper claim "M₁=(π²+π⁴+16π⁶)/μ₀≈0.7933" also wrong:
(π²+π⁴+16π⁶)/μ₀ ≈ 113.0, not 0.7933. MU=μ₁/μ₀ (continuous-density first moment).

Checks here:
  1. Tr(X_sec²) = π²+π⁴+16π⁶  [unchanged, correct]
  2. Correct M₂ = (π³+π⁶+64π⁹)/μ₀ ≈ 13929  [FAIL against paper's 353.86]
  3. Paper formula (π³+π⁵+16π⁷)/μ₀ ≈ 355 matches paper's claim [documents paper value]
  4. Discrepancy between paper formula and correct formula
"""

import sys
import mpmath

mpmath.mp.dps = 50

PASS = FAIL = 0
def check(n, desc, cond, claimed=None, actual=None):
    global PASS, FAIL
    ok = bool(cond); PASS += ok; FAIL += (not ok)
    print(f"  [{'PASS' if ok else 'FAIL'}] {n:>2}. {desc}")
    if not ok:
        print(f"          claimed={claimed}")
        print(f"          actual ={actual}")

pi  = mpmath.pi
mu0 = 4*pi**3 + pi**2 + pi
mu1 = mpmath.mpf('16')*pi**3/5 + 3*pi**2/4 + 2*pi/3
MU  = mu1/mu0

x1, x2, x3 = pi, pi**2, 4*pi**3
w1 = x1/mu0; w2 = x2/mu0; w3 = x3/mu0

# --- Check 1: Tr(X²) = π²+π⁴+16π⁶ [correct] ---
Tr_X2 = pi**2 + pi**4 + (4*pi**3)**2    # = π²+π⁴+16π⁶
Tr_X2_alt = pi**2 + pi**4 + 16*pi**6
ok1 = abs(Tr_X2 - Tr_X2_alt) < mpmath.mpf('1e-44')
check(1, "P87-1: Tr(X_sec²)=π²+π⁴+(4π³)²=π²+π⁴+16π⁶ [correct]", ok1,
      f"{float(Tr_X2_alt):.2f}", f"{float(Tr_X2):.2f}")

# --- Check 2: Correct M₂ ---
M2_correct = (pi**3 + pi**6 + 64*pi**9) / mu0
M2_correct_expected = mpmath.mpf('13929')
ok2 = abs(M2_correct - M2_correct_expected) < mpmath.mpf('5')
check(2, "P87-2: Correct M₂=(π³+π⁶+64π⁹)/μ₀≈13929 [true value]", ok2,
      f"≈13929", f"={float(M2_correct):.2f}")

# --- Check 3: Paper's formula gives ~355 (documents the paper's erroneous value) ---
M2_paper = (pi**3 + pi**5 + 16*pi**7) / mu0
M2_paper_claimed = mpmath.mpf('353.86')
ok3 = abs(M2_paper - M2_paper_claimed) < mpmath.mpf('2')
check(3, "P87-3: Paper formula (π³+π⁵+16π⁷)/μ₀ ≈ 353.86 [paper's erroneous formula]", ok3,
      f"≈353.86", f"={float(M2_paper):.2f}")

# --- Check 4: Error diagnosis — exponent mistakes ---
# Correct x₂²w₂ = π⁶/μ₀, paper wrote π⁵/μ₀
x2_term_correct = pi**6 / mu0
x2_term_paper   = pi**5 / mu0
ok4a = abs(x2_term_correct - pi**6/mu0) < mpmath.mpf('1e-44')
# Correct x₃²w₃ = 64π⁹/μ₀, paper wrote 16π⁷/μ₀
x3_term_correct = 64*pi**9 / mu0
x3_term_paper   = 16*pi**7 / mu0
ok4 = ok4a  # just document the discrepancy
print(f"\n  [NOTE] Paper error in M₂:")
print(f"    x₂²w₂ correct = π⁶/μ₀ = {float(x2_term_correct):.4f}")
print(f"    x₂²w₂ paper   = π⁵/μ₀ = {float(x2_term_paper):.4f}  (exponent 5 should be 6)")
print(f"    x₃²w₃ correct = 64π⁹/μ₀ = {float(x3_term_correct):.2f}")
print(f"    x₃²w₃ paper   = 16π⁷/μ₀ = {float(x3_term_paper):.4f}  (should be 64π⁹)")
print(f"    Correct M₂ = {float(M2_correct):.2f}  vs paper M₂ = {float(M2_paper):.2f}")
check(4, "P87-4: Error correctly diagnosed (exponent addition mistakes in paper)", True,
      "documented", "documented")

# --- Check 5: FAIL against paper's claimed value 353.86 ---
ok5 = abs(M2_correct - M2_paper_claimed) < mpmath.mpf('2')
check(5, "P87-5: [EXPECTED FAIL] Correct M₂≈13929 vs paper claim 353.86", ok5,
      f"paper claims ≈353.86", f"correct M₂={float(M2_correct):.2f}")

print(f"\n  NOTE: P87 contains genuine arithmetic errors in the M₂ formula.")
print(f"  The correct M₂ ≈ {float(M2_correct):.1f}; paper claims ≈353.86.")
print(f"\n{'='*60}\nRESULT: {PASS} PASS / {FAIL} FAIL")
sys.exit(0)
