"""
verify_P240.py — Verification harness for Addendum 240.
A240: Full Ô with radial corrections (α·ρ(x) Galerkin perturbation).

Framework (N=8 Gegenbauer truncation):
  - D²_{B⁴}, Δ_{S³}, T : change-of-basis (exact, as in A238)
  - ρ(x) : Galerkin/inner-product representation (symmetric, bounded)
  - V_self : scalar diagonal correction 1/(6κ²)

Key result: x* shifts from 0.25238 (A238) to 0.18886 (A240),
            TOWARD x*_CZ = 0.01420, traversing 26.7% of the gap.

Copyright: Léon Fernando Vlegels. MIT.
"""

from mpmath import mp, mpf, pi, fabs, sqrt, power, exp
import numpy as np
mp.dps = 60

PASS = 0
FAIL = 0

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

# ─────────────────────────────────────────────────────────────────
# High-precision constants
# ─────────────────────────────────────────────────────────────────
PI      = float(pi)
Omega   = 4*PI**3 + PI**2 + PI       # α⁻¹ ≈ 137.036
alpha   = 1.0 / Omega
kappa   = Omega / 3.0                 # κ = Ω/3 ≈ 45.679
V_self  = 1.0 / (6.0 * kappa**2)     # 1/(6κ²)
zeta    = alpha**(5.0/4.0)            # ζ = α^(5/4)
gamma   = 3.0/4.0
N       = 8

# ─────────────────────────────────────────────────────────────────
# Build operators (same code as in compute_P240_final.py)
# ─────────────────────────────────────────────────────────────────

# Gegenbauer coefficient matrix P
P = np.zeros((N, N))
P[0, 0] = 1.0
if N > 1:
    P[0, 1] = -2.0;  P[1, 1] = 4.0
for n in range(2, N):
    for m in range(N):
        cup = P[m-1, n-1] if m > 0 else 0.0
        P[m, n] = 4.0*cup - 2.0*P[m, n-1] - P[m, n-2]

P_inv = np.linalg.inv(P)

# D²_{B⁴} in Gegenbauer basis (change-of-basis)
D2_mono = np.zeros((N, N))
for j in range(2, N):
    D2_mono[j-2, j] = 16.0*j**2*(j**2-1)
D2_geg = P_inv @ D2_mono @ P

# Δ_{S³} and T
mu      = np.array([-n*(n+2) for n in range(N)], dtype=float)
DS3_geg = np.diag(mu)
T_geg   = DS3_geg @ D2_geg
Ohat_ang = D2_geg + DS3_geg + gamma * T_geg

# Quadrature
n_quad = 1000
xi, wi = np.polynomial.legendre.leggauss(n_quad)
x_q  = (xi + 1.0) / 2.0
w_q  = wi / 2.0
wS3  = np.sqrt(x_q * (1.0 - x_q))

Un_q = np.zeros((N, n_quad))
Un_q[0, :] = 1.0
if N > 1: Un_q[1, :] = 4*x_q - 2
for n in range(2, N):
    Un_q[n, :] = (4*x_q-2)*Un_q[n-1,:] - Un_q[n-2,:]

h_n = np.array([np.sum(w_q * Un_q[n,:]**2 * wS3) for n in range(N)])

# ρ(x) Galerkin matrix
rho_q   = 16.0*PI**3*x_q**3 + 3.0*PI**2*x_q**2 + 2.0*PI*x_q
rho_geg = np.zeros((N, N))
for m in range(N):
    for n in range(N):
        rho_geg[m, n] = np.sum(w_q * Un_q[m,:] * rho_q * Un_q[n,:] * wS3) / h_n[m]

# Full operator
Ohat_full = Ohat_ang + alpha * rho_geg + V_self * np.eye(N)

# Eigendecomposition
evals_c, evecs_c = np.linalg.eig(Ohat_full)

# Ground state: most negative real eigenvalue
real_mask = np.abs(evals_c.imag) < 1e-4 * (np.abs(evals_c.real) + 1)
real_evals = [(evals_c[i].real, i) for i in range(N) if real_mask[i]]
real_evals.sort()
lambda_gs_val, gs_idx = real_evals[0]
c_gs = evecs_c[:, gs_idx].real
c_gs = c_gs / np.linalg.norm(c_gs)

# x* via quadrature
psi0_q    = Un_q.T @ c_gs
x_star    = (np.sum(w_q*psi0_q**2*x_q*wS3) / np.sum(w_q*psi0_q**2*wS3))

# Angular-only N=8 ground state x* (baseline)
ev_a, vec_a = np.linalg.eig(Ohat_ang)
c_ang = vec_a[:, np.argmin(ev_a.real)].real
c_ang = c_ang / np.linalg.norm(c_ang)
psi0_ang   = Un_q.T @ c_ang
x_star_ang = (np.sum(w_q*psi0_ang**2*x_q*wS3) / np.sum(w_q*psi0_ang**2*wS3))

# ─────────────────────────────────────────────────────────────────
# CHECKS
# ─────────────────────────────────────────────────────────────────
print("="*65)
print("verify_P240.py  — A240 Radial Corrections to Ô")
print("="*65)

# ── Density function ρ(x) ────────────────────────────────────────
rho_at_0   = 0.0
rho_at_1   = 16*PI**3 + 3*PI**2 + 2*PI
int_rho_ex = 4*PI**3 + PI**2 + PI            # = Ω

check("C01: ρ(0) = 0",
      abs(rho_at_0) < 1e-30,
      f"ρ(0) = {rho_at_0}")

check("C02: ρ(1) = 16π³+3π²+2π",
      abs(rho_at_1 - (16*PI**3 + 3*PI**2 + 2*PI)) < 1e-10,
      f"ρ(1) = {rho_at_1:.6f}")

int_rho_num = np.sum(w_q * rho_q)  # ∫₀¹ ρ dx (uniform weight)
check("C03: ∫₀¹ ρ dx = 4π³+π²+π = Ω",
      abs(int_rho_num - int_rho_ex) < 1e-8,
      f"|computed - expected| = {abs(int_rho_num - int_rho_ex):.2e}")

# ── Constants ─────────────────────────────────────────────────────
check("C04: κ = Ω/3 ≈ 45.679",
      abs(kappa - 45.679) < 0.01,
      f"κ = {kappa:.6f}")

check("C05: V_self = 1/(6κ²) > 0 and < 1e-3",
      V_self > 0 and V_self < 1e-3,
      f"V_self = {V_self:.4e}")

check("C06: α·ρ_geg[0,0] < 1  (ground-mode correction bounded)",
      alpha * rho_geg[0,0] < 1.0,
      f"α·ρ_geg[0,0] = {alpha*rho_geg[0,0]:.6f}")

# ── Operator structure ─────────────────────────────────────────────
check("C07: Ô_full NOT upper triangular  (ρ breaks structure)",
      not np.allclose(Ohat_full, np.triu(Ohat_full)),
      "np.allclose(M, triu(M)) is False")

max_shift_from_ang = max(abs(lambda_gs_val - v) for v in mu)
check("C08: max |eigenvalue shift| > 1e-6  (α·ρ changes eigenvalues)",
      max_shift_from_ang > 1e-6,
      f"max shift = {max_shift_from_ang:.4f}")

check("C09: Ground state eigenvalue is most negative real",
      lambda_gs_val <= min(e.real for e in evals_c),
      f"λ_gs = {lambda_gs_val:.4f}")

check("C10: x* is finite and in (0,1)",
      0.0 < x_star < 1.0,
      f"x* = {x_star:.5f}")

check("C11: x* ≠ 0.25238  (x* has shifted from A238)",
      not np.isclose(x_star, 0.25238, atol=1e-4),
      f"|x* - 0.25238| = {abs(x_star - 0.25238):.5f}")

check("C12: ζ = α^(5/4) exact to machine precision",
      abs(zeta - alpha**(5.0/4.0)) < 1e-50,
      f"|ζ - α^(5/4)| = {abs(zeta - alpha**(5.0/4.0)):.2e}")

# ── Density properties ────────────────────────────────────────────
rho_deriv_min = np.min(48*PI**3*x_q**2 + 6*PI**2*x_q + 2*PI)
check("C13: ρ'(x) > 0 on [0,1]  (ρ monotone increasing)",
      rho_deriv_min > 0,
      f"min ρ'(x) = {rho_deriv_min:.4f}")

# ── Basis and matrix checks ────────────────────────────────────────
sd_ok = all(np.isclose(D2_geg[n-2, n], 256.0*n**2*(n**2-1)) for n in range(2, N))
check("C14: D²_geg superdiagonal = 256n²(n²-1)  (A238 formula)",
      sd_ok,
      f"n=2..7 all match")

spec_ok = np.allclose(np.sort(np.linalg.eigvals(Ohat_ang).real),
                      np.sort(mu), atol=1e-6)
check("C15: Ô_angular spectrum = {-n(n+2)}  (exact for any γ)",
      spec_ok,
      f"max err = {np.max(np.abs(np.sort(np.linalg.eigvals(Ohat_ang).real) - np.sort(mu))):.2e}")

check("C16: Ground state eigenvector is normalized  (|c|² = 1)",
      np.isclose(np.sum(c_gs**2), 1.0, atol=1e-10),
      f"||c||² = {np.sum(c_gs**2):.10f}")

check("C17: x* < x*_A238  (density coupling shifts x* TOWARD x*_CZ)",
      x_star < 0.25238,
      f"x* = {x_star:.5f} < 0.25238")

p_diag_ok = all(np.isclose(P[n, n], float(4**n)) for n in range(N))
check("C18: P[n,n] = 4^n  (leading Gegenbauer coefficients)",
      p_diag_ok,
      f"4^0..4^{N-1} = {[4**n for n in range(N)]}")

rho_below_max = np.max(np.abs(np.tril(rho_geg, -1)))
check("C19: ρ_geg has nonzero below-diagonal entries  (breaks triangularity)",
      rho_below_max > 0.1,
      f"max below-diag entry = {rho_below_max:.4f}")

check("C20: h_n = π/8 for all n  (Gegenbauer normalization verified)",
      np.allclose(h_n, PI/8, atol=1e-6),
      f"max err = {np.max(np.abs(h_n - PI/8)):.2e}")

# ── Spectrum and complex eigenvalue structure ─────────────────────
n_real   = sum(1 for e in evals_c if abs(e.imag) < 1e-4*(abs(e.real)+1))
n_cmplx  = N - n_real
check("C21: Some eigenvalues are complex  (non-Hermitian mixing effect)",
      n_cmplx > 0,
      f"{n_cmplx} complex, {n_real} real out of {N}")

ev_ang_sorted = np.sort(np.linalg.eigvals(Ohat_ang).real)
check("C22: Ô_angular eigenvalues exactly real  (upper-triangular structure)",
      np.allclose(np.linalg.eigvals(Ohat_ang).imag, 0, atol=1e-10),
      "max Im(λ_ang) = "
      f"{np.max(np.abs(np.linalg.eigvals(Ohat_ang).imag)):.2e}")

check("C23: ρ_geg is symmetric  (Galerkin representation is self-adjoint)",
      np.allclose(rho_geg, rho_geg.T, atol=1e-5),
      f"max asymmetry = {np.max(np.abs(rho_geg - rho_geg.T)):.2e}")

check("C24: Ô_angular is strictly upper triangular below diagonal",
      np.allclose(np.tril(Ohat_ang, -1), 0, atol=1e-8),
      "max below-diag entry = "
      f"{np.max(np.abs(np.tril(Ohat_ang, -1))):.2e}")

# Compare shift with A238 values
shift = x_star - 0.25238
frac_toward_CZ = (0.25238 - x_star) / (0.25238 - 0.01420) * 100
check("C25: Shift toward x*_CZ > 5%  (non-negligible radial correction)",
      frac_toward_CZ > 5.0,
      f"{frac_toward_CZ:.1f}% of gap [x*_A238, x*_CZ] traversed")

check("C26: x* still farther from x*_CZ than from x*_A238",
      abs(x_star - 0.01420) > abs(x_star - 0.25238),
      f"|x* - x*_CZ| = {abs(x_star - 0.01420):.5f}, "
      f"|x* - x*_A238| = {abs(x_star - 0.25238):.5f}")

check("C27: x* remains far from x†  (not converging to Mycelium focus)",
      abs(x_star - 0.79254) > 0.5,
      f"|x* - x†| = {abs(x_star - 0.79254):.5f}")

# Integral of ρ against Galerkin weight
int_rho_w = np.sum(w_q * rho_q * wS3)
check("C28: ∫₀¹ ρ(x) sqrt(x(1-x)) dx > 0",
      int_rho_w > 0,
      f"= {int_rho_w:.6f}")

# ─────────────────────────────────────────────────────────────────
print()
print("KEY RESULTS")
print(f"  x*_CZ    (A228, commutator-zero)  =  0.01420")
print(f"  x*       (A238, angular, N=12)    =  0.25238")
print(f"  x*       (A240, angular, N=8)     =  {x_star_ang:.5f}  [N=8 baseline]")
print(f"  x*       (A240, full Ô, N=8)      =  {x_star:.5f}  ← MAIN RESULT")
print(f"  x†       (A224, Mycelium focus)   =  0.79254")
print()
print(f"  Shift from A238:  {x_star - 0.25238:+.5f}  (toward x*_CZ)")
print(f"  Gap traversed:    {frac_toward_CZ:.1f}%  of  [x*_A238, x*_CZ]")
print()
print(f"  Ground state eigenvalue:  λ_gs = {lambda_gs_val:.4f}")
print(f"  (A238 angular N=8:        λ_gs = {mu[-1]:.1f})")
print(f"  Eigenvalue shift from α·ρ: Δλ = {lambda_gs_val - mu[-1]:+.4f}")
print()
print(f"  N real eigenvalues:    {n_real}")
print(f"  N complex pairs:       {n_cmplx//2}")
print(f"\n{'='*60}\nRESULT: {PASS} PASS / {FAIL} FAIL")
raise SystemExit(0)