#!/usr/bin/env python3
"""
verify_P238.py — Addendum 238: Ô_angular eigenspectrum and eigentrajectory x*
N=12 Gegenbauer basis {U_0,...,U_11}. numpy-only (scipy not required).
Copyright: Léon Fernando Vlegels. License: MIT.
"""
from mpmath import mp, mpf, pi as mp_pi, fabs
import numpy as np
import sys

mp.dps = 60

PASS = 0; FAIL = 0; _N = 0

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

# ── Constants ──────────────────────────────────────────────────────────────────
N = 12
alpha_inv = 4*np.pi**3 + np.pi**2 + np.pi
alpha     = 1.0 / alpha_inv
gamma_    = 3.0 / 4.0

# Reference values from prior addenda
x_CZ     = (np.pi - 1.0) / (48.0 * np.pi)   # ≈ 0.014159…
x_dagger = 0.79254                            # Mycelium focus node (A224)

print(f"\n=== A238: Ô_angular Eigenspectrum  (N={N} Gegenbauer modes) ===\n")
print(f"  α⁻¹  = {alpha_inv:.8f}")
print(f"  γ    = {gamma_}")
print(f"  x*_CZ = {x_CZ:.8f}  (commutator-zero proxy, A228)")
print(f"  x†    = {x_dagger:.8f}  (Mycelium focus node, A224)")

# ── Build U_n(x) as monomial coefficient vectors ───────────────────────────────
# Recurrence: U_n(x) = (4x-2)*U_{n-1}(x) - U_{n-2}(x)
# coeffs[n, m] = coefficient of x^m in U_n(x)
coeffs = np.zeros((N, N), dtype=np.float64)
coeffs[0, 0] = 1.0                          # U_0 = 1
if N > 1:
    coeffs[1, 0] = -2.0                     # U_1 = 4x - 2
    coeffs[1, 1] =  4.0
for n in range(2, N):
    for m in range(N):
        v = -2.0 * coeffs[n-1, m]
        if m >= 1:
            v += 4.0 * coeffs[n-1, m-1]
        v -= coeffs[n-2, m]
        coeffs[n, m] = v

# Spot-check recurrence against known values
check("U_2 const term", abs(coeffs[2,0]-3)  < 1e-12)
check("U_2 x term", abs(coeffs[2,1]+16) < 1e-12)
check("U_2 x^2 term", abs(coeffs[2,2]-16) < 1e-12)
check("U_3 const term", abs(coeffs[3,0]+4)  < 1e-12)

# P[m, n] = coefficient of x^m in U_n(x)  →  P = coeffs.T
P = coeffs.T                                # shape (N, N), upper triangular

# ── D²_{B⁴} in monomial basis ─────────────────────────────────────────────────
# D²_mono[j-2, j] = 16*j²*(j²-1) for j >= 2, all other entries 0
D2_mono = np.zeros((N, N), dtype=np.float64)
for j in range(2, N):
    D2_mono[j-2, j] = 16.0 * j**2 * (j**2 - 1)

# ── D²_{B⁴} in Gegenbauer basis: [D²]_geg = P⁻¹ D²_mono P ───────────────────
P_inv  = np.linalg.inv(P)
D2_geg = P_inv @ D2_mono @ P

# ── Δ_{S³}: diagonal with μ_n = -n(n+2) ─────────────────────────────────────
mu       = np.array([-n*(n+2) for n in range(N)], dtype=np.float64)
Delta_S3 = np.diag(mu)

# ── T = Δ_{S³} ∘ D²_{B⁴} in Gegenbauer basis ────────────────────────────────
T = Delta_S3 @ D2_geg          # T[m,n] = μ_m * D2_geg[m,n]

# ── Ô_angular = D²_{B⁴} + Δ_{S³} + γ T ──────────────────────────────────────
O_ang = D2_geg + Delta_S3 + gamma_ * T

print(f"\n--- Matrix construction ---")
print(f"  D2_geg[0,2] = {D2_geg[0,2]:.4f}   (expect 3072 = 256·4·3 from superdiag)")
print(f"  T[1,3]      = {T[1,3]:.4f}   (expect -55296 = -3·256·9·8)")
diag_vals = ["%.1f" % O_ang[n,n] for n in range(6)]
print(f"  O_ang diagonal: {diag_vals}")

# ──────────────────────────────────────────────────────────────────────────────
# CHECK C01 — monomial D²_mono[0,2] = 192  (D²(x²) = 16·4·3·x⁰ = 192)
check("C01: D2_mono[0,2] = 192  (D²(x²) = 16·2²·(4-1)·x⁰)",
      abs(D2_mono[0, 2] - 192.0) < 1e-10)

# CHECK C02 — Δ_{S³} eigenvalues
check("C02a: μ_0 = 0",    abs(mu[0]) < 1e-12)
check("C02b: μ_1 = -3",   abs(mu[1] + 3.0) < 1e-12)
check("C02c: μ_2 = -8",   abs(mu[2] + 8.0) < 1e-12)
check("C02d: μ_3 = -15",  abs(mu[3] + 15.0) < 1e-12)

# CHECK C03 — T superdiagonal T[1,3] = μ_1 · D²_geg[1,3]
# D²_geg[1,3] = 256·3²·(3²-1) = 256·9·8 = 18432  (superdiag formula n=3)
T13_exp = -3.0 * 256.0 * 9.0 * 8.0     # = -55296
check(f"C03: T[1,3] = {T13_exp:.0f}",   abs(T[1, 3] - T13_exp) < 1e-4)

# CHECK C04 — Ô_angular is not symmetric (T term breaks symmetry)
check("C04: Ô_angular not symmetric",   not np.allclose(O_ang, O_ang.T, atol=1e-6))

# ── Diagonalize Ô_angular ─────────────────────────────────────────────────────
eigenvalues_c, eigenvectors_c = np.linalg.eig(O_ang)
eigenvalues  = np.real(eigenvalues_c)
eigenvectors = np.real(eigenvectors_c)

# Sort ascending
idx          = np.argsort(eigenvalues)
eigenvalues  = eigenvalues[idx]
eigenvectors = eigenvectors[:, idx]

print(f"\n--- Eigenspectrum (ascending) ---")
print(f"  {'n':>3}  {'λ_n':>14}  {'μ_n(diagonal)':>16}")
for i in range(N):
    print(f"  {i:>3}  {eigenvalues[i]:>14.6f}  {mu[N-1-i]:>16.6f}")

# Ground state (most negative eigenvalue)
psi0    = eigenvectors[:, 0]
lambda0 = eigenvalues[0]

# CHECK C05 — ground eigenvalue < 0
check(f"C05: Ground eigenvalue {lambda0:.4f} < 0",  lambda0 < 0.0)

# ── x* = <ψ₀|x|ψ₀>_w / <ψ₀|ψ₀>_w ─────────────────────────────────────────
# Use Gauss-Legendre quadrature on [0,1] with w(x) = sqrt(x(1-x))
nq   = 300
# Gauss-Legendre nodes/weights on [-1,1], mapped to [0,1]
xi, wi = np.polynomial.legendre.leggauss(nq)
xq   = 0.5 * (xi + 1.0)
wq   = 0.5 * wi

# Evaluate U_n at quadrature nodes via recurrence
Uq = np.zeros((N, nq))
Uq[0] = 1.0
if N > 1:
    Uq[1] = 4.0*xq - 2.0
for n in range(2, N):
    Uq[n] = (4.0*xq - 2.0) * Uq[n-1] - Uq[n-2]

wx = np.sqrt(xq * (1.0 - xq))          # weight function w(x) = √(x(1-x))

# Verify inner product normalization: <U_n, U_n>_w = π/4
norms_sq = np.sum(Uq**2 * wx * wq, axis=1)
print(f"\n--- Inner product check  <U_n,U_n>_w (expect π/4 = {np.pi/4:.8f}) ---")
for n in range(4):
    print(f"  n={n}: {norms_sq[n]:.8f}")

# x-matrix: X_raw[m,n] = ∫₀¹ U_m(x)·x·U_n(x)·w(x) dx
# ψ₀(x) = Σ_n c_n U_n(x); x* = (ψ₀, x·ψ₀)_w / (ψ₀, ψ₀)_w
# With eigenvector c = psi0, and (ψ₀,ψ₀)_w = Σ_n c_n² · π/4:
psi0_xq  = Uq.T @ psi0                              # ψ₀ at quadrature nodes
numerator   = np.sum(psi0_xq**2 * xq * wx * wq)    # <ψ₀|x|ψ₀>_w
denominator = np.sum(psi0_xq**2 * wx * wq)          # <ψ₀|ψ₀>_w
x_star      = numerator / denominator

print(f"\n--- Eigentrajectory ---")
print(f"  Ground eigenvalue λ₀     = {lambda0:.8f}")
print(f"  x* = <ψ₀|x|ψ₀>/<ψ₀|ψ₀>  = {x_star:.8f}")
print(f"  x*_CZ (A228)             = {x_CZ:.8f}  |Δ| = {abs(x_star - x_CZ):.6f}")
print(f"  x†   (A224)              = {x_dagger:.8f}  |Δ| = {abs(x_star - x_dagger):.6f}")

# Gegenbauer mode amplitudes in ground state
amp = psi0**2 / np.sum(psi0**2)
print(f"\n--- Ground state |c_n|² amplitudes ---")
for n in range(N):
    bar = "█" * int(amp[n]*50 + 0.5)
    print(f"  n={n:>2}: {amp[n]:.6f}  {bar}")

# CHECK C06 — x* ≠ 0
check("C06: x* ≠ 0",  abs(x_star) > 1e-8)

# CHECK C07 — x* within 0.5 of x*_CZ or x†
check("C07: |x* - x*_CZ| < 0.5  OR  |x* - x†| < 0.5",
      abs(x_star - x_CZ) < 0.5 or abs(x_star - x_dagger) < 0.5)

# CHECK C08 — exactly N eigenvalues
check(f"C08: len(eigenvalues) == {N}",  len(eigenvalues) == N)

# CHECK C09 — nilpotency: T^(N//2+1) ≈ 0
k_nil  = N // 2 + 1
T_pow  = np.linalg.matrix_power(T, k_nil)
max_T  = np.max(np.abs(T_pow))
print(f"\n  T^{k_nil} max-abs entry: {max_T:.3e}")
check(f"C09: T^{k_nil} ≈ 0  (nilpotency, A231/A236)",  max_T < 1e-3)

# CHECK C10 — superdiagonal of D²_geg: [D²]_geg[n-2,n] = 256·n²·(n²-1)
c10_ok = True
for n in range(2, N):
    exp_val = 256.0 * n**2 * (n**2 - 1)
    act_val = D2_geg[n-2, n]
    if abs(act_val - exp_val) > max(1.0, abs(exp_val)*1e-6):
        c10_ok = False
        print(f"  C10 fail n={n}: D2_geg[{n-2},{n}]={act_val:.4f}, expect {exp_val:.4f}")
check("C10: D²_geg superdiagonal = 256·n²·(n²-1)  (A231 Thm 2)", c10_ok)

# CHECK C11 — all eigenvalues finite and real
check("C11: All eigenvalues finite real",
      all(np.isfinite(eigenvalues)) and
      np.max(np.abs(np.imag(eigenvalues_c[idx]))) < 1e-8)

# CHECK C12 — diagonal of Ô_angular = μ_n (since D²_geg and T are strictly upper tri)
diag_ok = all(abs(O_ang[n,n] - mu[n]) < 1e-8 for n in range(N))
check("C12: Ô_angular diagonal = μ_n = −n(n+2)  (D²,T strictly upper-tri)", diag_ok)

# CHECK C13 — D2_geg strictly upper triangular (offset ≥ 2, per A231/A236)
strictly_ut = all(abs(D2_geg[m, n]) < 1e-10 for n in range(N) for m in range(n+1, N))
check("C13: D2_geg strictly upper triangular (offset ≥2)",  strictly_ut)

# CHECK C14 — μ_n consistency: eigenvalues of Ô_angular match diagonal entries
# (upper tri matrix has eigenvalues = diagonal)
mu_sorted = np.sort(mu)          # ascending
ev_rounded = np.round(eigenvalues, 4)
mu_rounded = np.round(mu_sorted, 4)
check("C14: Eigenvalues match diagonal {μ_n} (upper-tri property)",
      np.allclose(ev_rounded, mu_rounded, atol=1e-2))

# CHECK C15 — x* in (0,1) (it is a position expectation on [0,1])
check("C15: x* ∈ (0, 1)",  0.0 < x_star < 1.0)

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