"""
verify_P225.py — Verification suite for Addendum 225
"Peter-Weyl Spectral Decomposition of ρ on S³ — Boundary Mode Identification"

L. F. Vlegels, 22 May 2026
© Léon Fernando Vlegels. MIT License.

Assertions:
  P001–P006  Gegenbauer-Chebyshev basis: definition and eigenvalues
  P007–P013  Exact Gegenbauer coefficients of ρ
  P014–P018  a₃ = OMEGA_0 = π³/4 (central identity)
  P019–P024  Reconstruction error < machine precision
  P025–P030  Spectral energy distribution
  P031–P036  Near-miss: ρ(MU)·Var ≈ OMEGA_0 at 0.295%
  P037–P042  Boundary stratum conditional expectation ≠ α⁻¹
  P043–P046  Layer-sourcing: n=3 mode exclusively from bulk term

All high-precision assertions at mp.dps = 60.
Expected outcome: 46 assertions, all passing.
"""

import sys
import os
import math
import numpy as np

# ── mpmath ────────────────────────────────────────────────────────────────────
from mpmath import mp, mpf, pi as mpi, sqrt as msqrt, fabs, nstr, quad
mp.dps = 60

# ── Path setup ────────────────────────────────────────────────────────────────
REPO_ROOT = os.path.abspath(os.path.join(os.path.dirname(__file__), "..", "..", "..", ".."))
if REPO_ROOT not in sys.path:
    sys.path.insert(0, REPO_ROOT)

# ── Assertion infrastructure ──────────────────────────────────────────────────
PASS = 0
FAIL = 0

def check(name: str, condition: bool) -> None:
    global PASS, FAIL
    ok = bool(condition)
    PASS += ok
    FAIL += (not ok)
    print(f"  [{'PASS' if ok else 'FAIL'}] {PASS + FAIL:>2}. {name}")


# ── Constants (high precision) ────────────────────────────────────────────────
OMEGA    = 4*mpi**3 + mpi**2 + mpi       # ≈ 137.036
OMEGA_0  = mpi**3 / 4                    # kernel constant π³/4 ≈ 7.7516

FRAC_BULK = 4*mpi**3 / OMEGA
FRAC_BND  = mpi**2 / OMEGA
FRAC_EDGE = mpi / OMEGA

# Moments mₙ = ∫₀¹ xⁿ ρ(x) dx  (P02 Thm 3.3: 16π³/(n+4) + 3π²/(n+3) + 2π/(n+2))
def moment(n):
    n = mpf(n)
    return 16*mpi**3/(n+4) + 3*mpi**2/(n+3) + 2*mpi/(n+2)

M0 = moment(0)   # = Ω
M1 = moment(1)   # first moment
M2 = moment(2)   # second moment
MU = M1 / M0     # bulk mean  ≈ 0.7933

# Variance Var(x) = M2/M0 − MU²
VAR = M2/M0 - MU**2

# ─────────────────────────────────────────────────────────────────────────────
# P001–P006  Gegenbauer-Chebyshev basis: eigenvalues and orthogonality
# ─────────────────────────────────────────────────────────────────────────────
print("S1  P001–P006  Gegenbauer basis and S³ eigenvalues")

# S³ Laplacian eigenvalue for Chebyshev-U mode n: μ_n = -n(n+2)
mu = lambda n: -n*(n+2)

check("P001  μ_0 = 0  (constant mode is in the kernel)",
      mu(0) == 0)
check("P002  μ_1 = -3  (dipole mode)",
      mu(1) == -3)
check("P003  μ_2 = -8  (quadrupole mode)",
      mu(2) == -8)
check("P004  μ_3 = -15  (octupole mode)",
      mu(3) == -15)

# Eigenvalues are strictly decreasing
check("P005  μ_0 > μ_1 > μ_2 > μ_3  (strictly decreasing)",
      mu(0) > mu(1) > mu(2) > mu(3))

# The 4-term expansion terminates because ρ is degree 3
check("P006  ρ is degree 3 in x → expansion terminates at n=3 (exactly 4 terms)",
      True)  # structural assertion: polynomial of degree n → at most n+1 Gegenbauer terms

# ─────────────────────────────────────────────────────────────────────────────
# P007–P013  Exact Gegenbauer coefficients of ρ
# ─────────────────────────────────────────────────────────────────────────────
print("S2  P007–P013  Gegenbauer coefficients of ρ")

# In v = 2x-1 variable: ρ_v(v) = a₀U₀ + a₁U₁ + a₂U₂ + a₃U₃
# Coefficient matching from leading terms (A225 derivation):
#   v³ coefficient of ρ_v: 2π³   → U₃ leading = 8 → a₃ = 2π³/8 = π³/4
#   v² coefficient (after removing a₃ contribution): → a₂
#   etc.
#
# Alternatively, use inner product on [-1,1] with weight √(1-v²):
# aₙ = (2(n+1)/π) ∫₋₁¹ ρ_v(v) U_n(v) √(1-v²) dv

def rho_v(v):
    """ρ expressed in v=2x-1 variable."""
    x = (v + 1) / 2
    return 16*mpi**3 * x**3 + 3*mpi**2 * x**2 + 2*mpi * x

def chebU(n, v):
    """Chebyshev-U polynomial U_n(v)."""
    if n == 0:
        return mpf(1)
    elif n == 1:
        return 2*v
    elif n == 2:
        return 4*v**2 - 1
    elif n == 3:
        return 8*v**3 - 4*v
    else:
        raise ValueError(n)

def coeff_gegenbauer(n):
    """aₙ via inner product on S³ (zonal, weight √(1-v²))."""
    norm = mpf(n+1) * mpi / 2   # ∫₋₁¹ U_n² √(1-v²) dv = π/2 for U_n normalised?
    # Actually: ∫₋₁¹ U_n(v)² √(1-v²) dv = π/2  for all n ≥ 0
    # So aₙ = (2/π) ∫₋₁¹ ρ_v(v) U_n(v) √(1-v²) dv
    integrand = lambda v: rho_v(v) * chebU(n, v) * msqrt(1 - v**2)
    inner = quad(integrand, [-1, 1])
    return 2 * inner / mpi

a = [coeff_gegenbauer(n) for n in range(4)]

# P007–P010: coefficients are positive and ordered roughly as expected
check("P007  a₀ > 0  (scalar mode coefficient positive)",
      a[0] > 0)
check("P008  a₁ > 0  (dipole mode coefficient positive)",
      a[1] > 0)
check("P009  a₂ > 0  (quadrupole mode coefficient positive)",
      a[2] > 0)
check("P010  a₃ > 0  (octupole mode coefficient positive)",
      a[3] > 0)

# P011: a₀ is the largest coefficient (normalization contains the mean)
check("P011  a₀ > a₁ > a₂ > a₃  (decreasing with mode number)",
      a[0] > a[1] > a[2] > a[3])

# P012: a₃ is in the expected range [7.7, 7.8]
check("P012  a₃ ∈ [7.7, 7.8]  (≈ π³/4 ≈ 7.7516)",
      mpf('7.7') < a[3] < mpf('7.8'))

# P013: a₂ is in the expected range [48, 49]
check("P013  a₂ ∈ [48, 49]  (expected ≈ 48.36)",
      mpf('48') < a[2] < mpf('49'))

# ─────────────────────────────────────────────────────────────────────────────
# P014–P018  a₃ = OMEGA_0 = π³/4  (central identity)
# ─────────────────────────────────────────────────────────────────────────────
print("S3  P014–P018  Central identity: a₃ = OMEGA_0 = π³/4")

# Algebraic derivation: the leading v³ coefficient of ρ_v is
# 16π³·(1/2)³ = 16π³/8 = 2π³.
# U₃(v) has leading term 8v³.  So a₃ = 2π³/8 = π³/4.

a3_algebraic = mpi**3 / 4   # = OMEGA_0

check("P014  a₃ = π³/4 = OMEGA_0  (algebraic: leading v³ coefficient of ρ_v is 2π³; U₃ leading = 8)",
      fabs(a[3] - a3_algebraic) < mpf('1e-50'))

check("P015  OMEGA_0 = π³/4  (kernel constant definition)",
      fabs(OMEGA_0 - mpi**3/4) < mpf('1e-55'))

check("P016  a₃ = OMEGA_0  (inner-product ≡ algebraic, to 50 decimal places)",
      fabs(a[3] - OMEGA_0) < mpf('1e-50'))

# Numerical float check
PI = float(mpi)
a3_float = float(a[3])
OMEGA_0_float = PI**3 / 4
check("P017  a₃ = OMEGA_0 in float64 (residual < 1e-12)",
      abs(a3_float - OMEGA_0_float) < 1e-12)

# Cross-check: a₃ computed from algebraic formula vs quadrature
a3_quad = float(coeff_gegenbauer(3))
check("P018  Quadrature a₃ matches π³/4 to 1e-12 (numerical confirmation)",
      abs(a3_quad - OMEGA_0_float) < 1e-12)

# ─────────────────────────────────────────────────────────────────────────────
# P019–P024  Reconstruction error < machine precision
# ─────────────────────────────────────────────────────────────────────────────
print("S4  P019–P024  Reconstruction accuracy")

def rho_reconstruct(x_val):
    """Reconstruct ρ from Gegenbauer expansion at x ∈ [0,1]."""
    v = 2*x_val - 1
    return (a[0]*chebU(0, v) + a[1]*chebU(1, v) +
            a[2]*chebU(2, v) + a[3]*chebU(3, v))

def rho_exact(x_val):
    return 16*mpi**3 * x_val**3 + 3*mpi**2 * x_val**2 + 2*mpi * x_val

test_points = [mpf('0.0'), mpf('0.25'), mpf('0.5'), mpf('0.75'), mpf('1.0')]

for i, xp in enumerate(test_points):
    diff = fabs(rho_reconstruct(xp) - rho_exact(xp))
    check(f"P0{19+i}  Reconstruction at x={float(xp):.2f}: residual < 1e-45",
          diff < mpf('1e-45'))

# ─────────────────────────────────────────────────────────────────────────────
# P025–P030  Spectral energy distribution
# ─────────────────────────────────────────────────────────────────────────────
print("S5  P025–P030  Spectral energy ‖aₙ‖² distribution")

# Parseval (with Chebyshev-U norm π/2): total energy = Σ aₙ² · π/2
# Fractional energy of mode n: aₙ² / Σ aₙ²
energies = [float(a[n]**2) for n in range(4)]
total_energy = sum(energies)
fracs = [e / total_energy for e in energies]

check("P025  Spectral energy n=0: 40-55% of total",
      0.40 < fracs[0] < 0.55)
check("P026  Spectral energy n=1: 40-50% of total",
      0.40 < fracs[1] < 0.50)
check("P027  Spectral energy n=2: 5-10% of total",
      0.05 < fracs[2] < 0.10)
check("P028  Spectral energy n=3: 0.1-0.5% of total  (bulk mode is minor in energy)",
      0.001 < fracs[3] < 0.005)
check("P029  Modes n=0,1 together hold > 85% of energy",
      fracs[0] + fracs[1] > 0.85)
check("P030  Bulk mode n=3 (a₃=OMEGA_0) is spectrally subdominant despite OMEGA_0 exactness",
      fracs[3] < 0.01)

# ─────────────────────────────────────────────────────────────────────────────
# P031–P036  Near-miss: ρ(MU)·Var ≈ OMEGA_0 at 0.295%
# ─────────────────────────────────────────────────────────────────────────────
print("S6  P031–P036  Near-miss ρ(MU)·Var ≈ OMEGA_0 (from HUP research file)")

rho_MU = rho_exact(MU)    # ρ evaluated at the bulk mean
product = rho_MU * VAR    # ρ(MU) × Var(x)

check("P031  ρ(MU) > 0  (density is positive at bulk mean)",
      rho_MU > 0)
check("P032  Var(x) > 0  (variance is positive)",
      VAR > 0)
check("P033  ρ(MU)·Var is close to OMEGA_0 (within 1%)",
      fabs(product - OMEGA_0) / OMEGA_0 < mpf('0.01'))
check("P034  ρ(MU)·Var is NOT exactly OMEGA_0 (residual > 1e-6)",
      fabs(product - OMEGA_0) > mpf('1e-6'))

# Near-miss is at ~0.295%:
rel_diff = fabs(product - OMEGA_0) / OMEGA_0
check("P035  Near-miss relative error in [0.002, 0.004]  (≈ 0.295%)",
      mpf('0.002') < rel_diff < mpf('0.004'))

check("P036  ρ(MU)·Var > OMEGA_0  (near-miss overshoots)",
      product > OMEGA_0)

# ─────────────────────────────────────────────────────────────────────────────
# P037–P042  Boundary stratum conditional expectation ≠ α⁻¹
# ─────────────────────────────────────────────────────────────────────────────
print("S7  P037–P042  Boundary stratum expectation ≠ α⁻¹")

# Boundary stratum: x in [FRAC_BULK, 1] (the top (f_bnd + f_edge) of ρ-measure)
x_lo = FRAC_BULK   # ≈ 0.905

numerator   = quad(lambda x: rho_exact(x), [x_lo, mpf(1)])  # ∫_{x_lo}^1 ρ dx
denominator = M0  # total: Ω

bnd_exp = numerator   # this is ∫_{x_lo}^1 ρ dx (NOT normalised — it's the partial integral)
CODATA_ALPHA_INV = mpf('137.035999084')

check("P037  ∫_{FRAC_BULK}^1 ρ dx > 0  (boundary stratum has nonzero ρ-mass)",
      bnd_exp > 0)
check("P038  ∫_{FRAC_BULK}^1 ρ dx < Ω  (partial integral is less than total)",
      bnd_exp < M0)

# The partial integral ≠ α⁻¹:
check("P039  ∫_{FRAC_BULK}^1 ρ dx ≠ α⁻¹_CODATA  (boundary stratum integral does not equal α⁻¹)",
      fabs(bnd_exp - CODATA_ALPHA_INV) > mpf('1'))

# The partial integral is about 43.93, well above α⁻¹:
check("P040  ∫_{FRAC_BULK}^1 ρ dx in [40, 50]",
      mpf('40') < bnd_exp < mpf('50'))

# The x* threshold (where ∫_{x*}^1 ρ dx = α⁻¹) sits in edge stratum
# x* ≈ 0.0097 < FRAC_EDGE ≈ 0.023 — check numerically
# (simple bisection for the threshold)
def partial_integral_from(x_lo_val):
    return quad(lambda x: rho_exact(x), [x_lo_val, mpf(1)])

# Binary search for x* such that ∫_{x*}^1 ρ dx = α⁻¹
lo, hi = mpf('0'), mpf('0.1')
for _ in range(60):
    mid = (lo + hi) / 2
    if partial_integral_from(mid) > CODATA_ALPHA_INV:
        lo = mid
    else:
        hi = mid
x_star = (lo + hi) / 2

check("P041  x* (threshold where ∫_{x*}^1 ρ dx = α⁻¹) is in edge stratum [0, FRAC_EDGE]",
      mpf('0') < x_star < FRAC_EDGE)

# The gap Δ = Ω − α⁻¹
DELTA = OMEGA - CODATA_ALPHA_INV
check("P042  Gap Δ = Ω − α⁻¹ ∈ [3.0e-4, 3.1e-4]  (2.22 ppm)",
      mpf('3.0e-4') < DELTA < mpf('3.1e-4'))

# ─────────────────────────────────────────────────────────────────────────────
# P043–P046  Layer-sourcing: n=3 mode exclusively from bulk term
# ─────────────────────────────────────────────────────────────────────────────
print("S8  P043–P046  Layer sourcing of n=3 (bulk) mode")

# The n=3 Gegenbauer mode of ρ is sourced exclusively by the 16π³x³ term:
# If we decompose ρ = ρ_edge + ρ_bnd + ρ_bulk = 2πx + 3π²x² + 16π³x³,
# then a₃(ρ_edge) = a₃(ρ_bnd) = 0, a₃(ρ_bulk) = π³/4.

def a3_of_term(poly_func):
    """a₃ coefficient of an arbitrary polynomial term."""
    integrand = lambda v: poly_func(v) * chebU(3, v) * msqrt(1 - v**2)
    return 2 * quad(integrand, [-1, 1]) / mpi

# Edge term: 2πx = π(v+1) in v
a3_edge = a3_of_term(lambda v: mpi * (v + 1))
# Boundary term: 3π²x² = (3π²/4)(v+1)² in v
a3_bnd  = a3_of_term(lambda v: (3*mpi**2/4) * (v + 1)**2)
# Bulk term: 16π³x³ = 2π³(v+1)³ in v
a3_bulk = a3_of_term(lambda v: 2*mpi**3 * (v + 1)**3)

check("P043  n=3 coefficient of edge term 2πx = 0  (edge does not source bulk mode)",
      fabs(a3_edge) < mpf('1e-45'))
check("P044  n=3 coefficient of boundary term 3π²x² = 0  (boundary does not source bulk mode)",
      fabs(a3_bnd) < mpf('1e-45'))
check("P045  n=3 coefficient of bulk term 16π³x³ = π³/4 = OMEGA_0  (bulk exclusively sources n=3)",
      fabs(a3_bulk - OMEGA_0) < mpf('1e-45'))
check("P046  a₃(edge) + a₃(boundary) + a₃(bulk) = a₃(ρ) = OMEGA_0  (layer additivity)",
      fabs(a3_edge + a3_bnd + a3_bulk - a[3]) < mpf('1e-45'))

# ─────────────────────────────────────────────────────────────────────────────
# Summary
# ─────────────────────────────────────────────────────────────────────────────

print(f"\nKey values:")
print(f"  Ω       = {nstr(OMEGA, 15)}")
print(f"  OMEGA_0 = {nstr(OMEGA_0, 15)}")
print(f"  a₃      = {nstr(a[3], 15)}  (should equal OMEGA_0)")
print(f"  a₂      = {nstr(a[2], 10)}")
print(f"  a₁      = {nstr(a[1], 10)}")
print(f"  a₀      = {nstr(a[0], 10)}")
print(f"  MU      = {nstr(MU, 10)}")
print(f"  Var(x)  = {nstr(VAR, 10)}")
print(f"  ρ(MU)·Var = {nstr(rho_MU*VAR, 10)}  (vs OMEGA_0 = {nstr(OMEGA_0, 10)})")
print(f"  x*      = {nstr(x_star, 8)}  (threshold for α⁻¹ partial integral)")
print(f"  FRAC_EDGE = {nstr(FRAC_EDGE, 8)}")
print(f"  Spectral fracs: {[f'{f:.4f}' for f in fracs]}")

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