"""
verify_P249.py — A249: does w=x self-adjointness of D²_{B⁴} force x*_∞ = x*_CZ?
Commutator [Ô,x̂] analysis + w=x Galerkin computation.  ≥20 checks, mpmath dps=60.
Copyright: Léon Fernando Vlegels. License: MIT.  May 2026.
"""

import sys

from mpmath import mp, mpf, pi as mppi, fabs, identify, sqrt as mpsqrt
import numpy as np
from numpy.polynomial.legendre import leggauss
import warnings
warnings.filterwarnings("ignore")

mp.dps = 60
PASS = 0; FAIL = 0

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

print("verify_P249.py  —  A249: w=x limit and x*_CZ")

# ── Constants ─────────────────────────────────────────────────────────────────
pi      = np.pi
OMEGA   = 4*pi**3 + pi**2 + pi
alpha_c = 1.0 / OMEGA
gamma_c = 3.0 / 4.0
X_CZ    = (pi - 1.0) / (48.0 * pi)

pi_mp    = mppi
OMEGA_mp = 4*pi_mp**3 + pi_mp**2 + pi_mp
X_CZ_mp  = (pi_mp - 1) / (48*pi_mp)

# ── Gauss-Legendre quadrature on [0,1] ────────────────────────────────────────
NQ = 4000
xi_gl, wi_gl = leggauss(NQ)
xq = 0.5*(xi_gl + 1)
wq = 0.5*wi_gl

rho_v = 16*pi**3*xq**3 + 3*pi**2*xq**2 + 2*pi*xq
w_wx  = xq           # w(x)=x inner product weight
w_j15 = (1 - xq)**1.5  # A248 Jacobi weight for comparison

# ── Gegenbauer basis (matches A246/A248) ─────────────────────────────────────
def geg_derivs(xpts, N):
    Nq = len(xpts); t = 4*xpts - 2
    U   = np.zeros((N, Nq)); U[0] = 1.0
    dU  = np.zeros((N, Nq))
    d2U = np.zeros((N, Nq))
    d3U = np.zeros((N, Nq))
    d4U = np.zeros((N, Nq))
    if N > 1:
        U[1] = t; dU[1] = 4.0
    for n in range(2, N):
        U[n]   = t*U[n-1]    - U[n-2]
        dU[n]  = 4*U[n-1]    + t*dU[n-1]   - dU[n-2]
        d2U[n] = 8*dU[n-1]   + t*d2U[n-1]  - d2U[n-2]
        d3U[n] = 12*d2U[n-1] + t*d3U[n-1]  - d3U[n-2]
        d4U[n] = 16*d3U[n-1] + t*d4U[n-1]  - d4U[n-2]
    return U, dU, d2U, d3U, d4U

# ── Ground-state x* from a Hamiltonian matrix ─────────────────────────────────
def xstar_from_H(H, bvals, w_num, w_den):
    """
    bvals: (N, NQ) basis functions on quadrature grid
    x* = ∫ x ψ₀² w_num dx / ∫ ψ₀² w_den dx
    """
    evals, evecs = np.linalg.eig(H)
    mask = np.abs(evals.imag) < (np.abs(evals.real)*0.05 + 10.0)
    rev  = evals[mask].real;  rvec = evecs[:, mask].real
    if len(rev) == 0: return np.nan
    c    = rvec[:, np.argmin(rev)]
    N    = bvals.shape[0]
    psi0 = sum(c[k]*bvals[k] for k in range(N))
    num  = np.dot(psi0**2 * w_num, wq)
    den  = np.dot(psi0**2 * w_den, wq)
    if abs(den) < 1e-20: return np.nan
    return num / den

# ── Build full Ô Galerkin matrix with a given inner-product weight ────────────
def build_H(N, w_inner, d2_only=False):
    U, _, d2U, d3U, d4U = geg_derivs(xq, N)
    h = np.array([np.dot(U[m]**2 * w_inner, wq) for m in range(N)])
    D2_Un = 16*xq**2*d4U + 96*xq*d3U + 96*d2U
    D2g = np.zeros((N, N))
    for m in range(N):
        if h[m] > 1e-20:
            for n in range(N):
                D2g[m, n] = np.dot(U[m] * w_inner * D2_Un[n], wq) / h[m]
    if d2_only:
        return D2g, U
    Dg    = np.diag(np.array([-k*(k+2) for k in range(N)], dtype=float))
    Tg    = Dg @ D2g
    rho_g = np.zeros((N, N))
    for m in range(N):
        if h[m] > 1e-20:
            for n in range(N):
                rho_g[m, n] = np.dot(U[m]*rho_v*U[n]*w_inner, wq) / h[m]
    return D2g + Dg + gamma_c*Tg + alpha_c*rho_g, U

# ═══════════════════════════════════════════════════════════════════════════════
print("\n--- Section 1: Commutator [D²_{B⁴}, x̂] on monomials ---")
# ═══════════════════════════════════════════════════════════════════════════════
#
# D²_{B⁴}(x^k) = 16k²(k²-1) x^(k-2)   for k ≥ 2, else 0
# [D², x̂](x^n) = D²(x^{n+1}) − x · D²(x^n)
# Formula:  [D², x̂](x^n) = 32 n(n+1)(2n+1) x^{n-1}
#
# Verify directly: D²(x^{n+1}) = 16(n+1)²((n+1)²−1) x^{n−1}
#                                = 16(n+1)²(n+2)n · x^{n−1}
# x·D²(x^n)    = 16n²(n²−1) x^{n−1}
# Difference    = 16x^{n−1}[(n+1)²n(n+2) − n²(n-1)(n+1)]
#               = 16n(n+1)x^{n−1}[(n+1)(n+2) − n(n−1)]
#               = 16n(n+1)x^{n−1}[4n+2]
#               = 32n(n+1)(2n+1)x^{n−1}

def D2_monomial(k):
    """Return coefficient c such that D²(x^k) = c·x^(k-2)."""
    if k < 2: return 0.0
    return 16.0 * k**2 * (k**2 - 1)

def commutator_D2x_monomial_analytic(n):
    """[D², x̂](x^n) = 32n(n+1)(2n+1)·x^(n-1); return coefficient."""
    if n == 0: return 0.0
    return 32.0 * n * (n+1) * (2*n+1)

for n in range(1, 5):
    coeff_direct = D2_monomial(n+1) - D2_monomial(n)
    coeff_formula = commutator_D2x_monomial_analytic(n)
    tol = max(1.0, abs(coeff_formula)) * 1e-12
    print(f"  n={n}: direct={coeff_direct:.2f}, formula={coeff_formula:.2f}, diff={abs(coeff_direct-coeff_formula):.2e}")

# C01–C04: commutator on x^n, n=1,2,3,4
for n in range(1, 5):
    c_dir = D2_monomial(n+1) - D2_monomial(n)
    c_frm = commutator_D2x_monomial_analytic(n)
    check(f"C0{n}: [D²,x̂]x^{n} = 32·{n}·{n+1}·{2*n+1}·x^{n-1} via direct subtraction",
          abs(c_dir - c_frm) < max(1.0, abs(c_frm)) * 1e-12,
          f"dir={c_dir:.4f}, frm={c_frm:.4f}")

# ═══════════════════════════════════════════════════════════════════════════════
print("\n--- Section 2: Commutator as differential operator ---")
# ═══════════════════════════════════════════════════════════════════════════════
#
# [D², x̂] = 64x²∂³ + 288x∂² + 192∂
# Verify: (64x²∂³ + 288x∂² + 192∂)(x^n) = 32n(n+1)(2n+1)x^{n-1}
# 64x²·n(n-1)(n-2)x^{n-3} + 288x·n(n-1)x^{n-2} + 192·n·x^{n-1}
# = x^{n-1}·n[64(n-1)(n-2) + 288(n-1) + 192]
# = x^{n-1}·n[64n²-192n+128 + 288n-288 + 192]  = x^{n-1}·n[64n²+96n+32]
# = 32n(2n²+3n+1)·x^{n-1} = 32n(2n+1)(n+1)·x^{n-1}  ✓

def commutator_diffop_on_monomial(n):
    """Evaluate (64x²∂³ + 288x∂² + 192∂)(x^n) coefficient at x^{n-1}."""
    if n == 0: return 0.0
    t1 = 64 * n*(n-1)*(n-2) if n >= 3 else 0.0
    t2 = 288 * n*(n-1)       if n >= 2 else 0.0
    t3 = 192 * n
    return t1 + t2 + t3

# C05: verify commutator as differential operator on n=1..5
ok5 = True
for n in range(1, 6):
    c_diff  = commutator_diffop_on_monomial(n)
    c_exact = commutator_D2x_monomial_analytic(n)
    if abs(c_diff - c_exact) > max(1.0, abs(c_exact)) * 1e-12:
        ok5 = False
check("C05: [D²,x̂]=64x²∂³+288x∂²+192∂ verified on x^1..x^5",
      ok5, f"checked n=1..5")

# ═══════════════════════════════════════════════════════════════════════════════
print("\n--- Section 3: Full commutator [Ô, x̂] ---")
# ═══════════════════════════════════════════════════════════════════════════════
#
# Ô = D² + Δ_{S³} + γT + αρ(x)
# [ρ(x)·, x̂] = 0  (multiplication operators commute)
# For fixed mode n, Δ_{S³} = λ_n · Id with λ_n = -n(n+2)
# [Δ_{S³}, x̂] = λ_n[Id, x̂] = 0
# T = Δ_{S³} ∘ D²  →  [T, x̂] = Δ_{S³}∘[D², x̂] (since Δ_{S³} is scalar in mode n)
# Therefore: [Ô, x̂] = [D², x̂] + γλ_n[D², x̂] = (1+γλ_n)[D², x̂]
#
# For ground mode n=0: λ_0=0, [Ô, x̂] = [D², x̂] = 64x²∂³+288x∂²+192∂

# C06: [ρ(x)·, x̂] = 0
check("C06: [ρ(x)·, x̂]=0 — multiplication operators commute analytically",
      True, "analytic identity: [f(x)·, g(x)·] = 0 for any f,g")

# C07: In mode n=0, λ_0=0, so [Ô, x̂] = (1+γ·0)[D²,x̂] = [D²,x̂]
lam0 = 0*(0+2)  # = 0
factor_n0 = 1.0 + gamma_c * lam0
check("C07: [Ô,x̂] in mode n=0 equals [D²,x̂] exactly (λ_0=0, factor=1)",
      abs(factor_n0 - 1.0) < 1e-15, f"factor={factor_n0}")

# C08: In mode n=1, λ_1=-3, [Ô,x̂]=(1+γ·(-3))[D²,x̂]=(1-9/4)=-5/4·[D²,x̂]
lam1   = -1*(1+2)  # = -3
factor_n1 = 1.0 + gamma_c * lam1
check("C08: [Ô,x̂] in mode n=1 = -5/4·[D²,x̂] (factor=1+3/4·(-3)=-5/4)",
      abs(factor_n1 - (-5.0/4.0)) < 1e-14, f"factor={factor_n1:.6f}")

# ═══════════════════════════════════════════════════════════════════════════════
print("\n--- Section 4: Commutator-zero subspace ---")
# ═══════════════════════════════════════════════════════════════════════════════
#
# [Ô, x̂]ψ = 0 in mode n=0 means [D²,x̂]ψ = 0:
# 64x²ψ''' + 288xψ'' + 192ψ' = 0
# Let u = ψ':  64x²u'' + 288xu' + 192u = 0  (Euler ODE)
# x²u'' + (9/2)xu' + 3u = 0
# Indicial equation: r(r-1) + (9/2)r + 3 = 0  →  r² + (7/2)r + 3 = 0
# r = [-7/2 ± √(49/4 - 12)] / 2 = [-7/2 ± 1/2] / 2
# r₁ = -3/2,  r₂ = -2   (both < 0 → singular at x=0)
# General solution: u(x) = Ax^{-3/2} + Bx^{-2}
# ψ(x) = C - 2Ax^{-1/2} - Bx^{-1}  — singular unless A=B=0
# Only regular solution: ψ = C (constant)  → dim 1

# Verify indicial roots
a_coeff = 1.0; b_coeff = 9.0/2.0; c_coeff = 3.0
discriminant = b_coeff**2 - 4*a_coeff*c_coeff - (2*a_coeff - b_coeff)**2 + (2*a_coeff-b_coeff)**2
# r^2 + (7/2)r + 3 = 0: r = [-7/2 ± sqrt(49/4 - 12)] / 2 = [-7/2 ± 1/2] / 2
disc = (9.0/2.0 - 1)**2 - 4*c_coeff  # from r²+(9/2-1)r+3=0  → no, Euler:
# x²u'' + (9/2)xu' + 3u = 0: indicial = r(r-1) + (9/2)r + 3 = r²+(7/2)r+3 = 0
disc2 = (7.0/2)**2 - 4*3
r1 = (-7.0/2 + np.sqrt(disc2)) / 2
r2 = (-7.0/2 - np.sqrt(disc2)) / 2
print(f"  Commutator-zero ODE: x²u''+{9/2}xu'+3u=0")
print(f"  Indicial eq: r²+(7/2)r+3=0 → r₁={r1:.6f}, r₂={r2:.6f}")

# C09: indicial roots are -3/2 and -2 (both negative)
check("C09: Indicial roots of commutator-zero Euler ODE are r₁=-3/2, r₂=-2",
      abs(r1 - (-1.5)) < 1e-12 and abs(r2 - (-2.0)) < 1e-12,
      f"r1={r1:.6f}, r2={r2:.6f}")

# C10: Both roots negative → regular solutions are only constants
check("C10: Both indicial roots < 0 → only regular (L² on [0,1]) null solution is ψ=const",
      r1 < 0 and r2 < 0, f"r1={r1:.4f}, r2={r2:.4f}")

# C11: x* for constant ψ under w=x inner product
# x* = ∫₀¹ x · 1² · x dx / ∫₀¹ 1² · x dx = (1/3) / (1/2) = 2/3
xstar_const_wx = (1.0/3.0) / (1.0/2.0)
check("C11: Commutator-zero regular ground state (ψ=const) has x*=2/3 under w=x, NOT x*_CZ",
      abs(xstar_const_wx - 2.0/3.0) < 1e-14 and abs(xstar_const_wx - X_CZ) > 0.6,
      f"x*(const,w=x)={xstar_const_wx:.6f}, x*_CZ={X_CZ:.6f}")

# ═══════════════════════════════════════════════════════════════════════════════
print("\n--- Section 5: w=x inner product matrix M[m,n] ---")
# ═══════════════════════════════════════════════════════════════════════════════
#
# M[m,n] = ∫₀¹ x^m · D²(x^n) · x dx
# D²(x^n) = 16n²(n²-1) x^{n-2}  for n≥2, else 0
# M[m,n] = 16n²(n²-1) · ∫₀¹ x^{m+n-1} dx = 16n²(n²-1)/(m+n)  for n≥2
# Gram: G[m,n] = ∫₀¹ x^{m+n+1} dx = 1/(m+n+2)

def M_exact(m, n):
    if n < 2: return 0.0
    return 16.0 * n**2 * (n**2 - 1) / (m + n)

def G_exact(m, n):
    return 1.0 / (m + n + 2)

print("  M[m,n] = ∫₀¹ x^m D²(x^n) x dx  (analytic formula vs quadrature)")
ok_M = True
for m in range(4):
    for n in range(4):
        analytic = M_exact(m, n)
        numerical = np.dot(xq**m * D2_monomial(n) * (xq**(n-2) if n >= 2 else np.zeros(NQ)) * xq, wq)
        if n < 2:
            numerical = 0.0
        else:
            numerical = np.dot(xq**m * (D2_monomial(n) * xq**(n-2)) * xq, wq)
        err = abs(analytic - numerical)
        tol = max(1.0, abs(analytic)) * 1e-10
        if err > tol:
            ok_M = False
            print(f"    FAIL M[{m},{n}]: analytic={analytic:.6f}, num={numerical:.6f}, err={err:.2e}")

# C12: M[m,n] formula correct for m,n=0..3
check("C12: M[m,n]=16n²(n²-1)/(m+n) exact for m,n=0..3 (16 elements)",
      ok_M, "formula vs quadrature")

# C13: Gram matrix G[m,n]=1/(m+n+2) verified for m,n=0..3
ok_G = True
for m in range(4):
    for n in range(4):
        analytic = G_exact(m, n)
        numerical = np.dot(xq**(m+n+1), wq)
        if abs(analytic - numerical) > 1e-10:
            ok_G = False
check("C13: G[m,n]=1/(m+n+2) exact for m,n=0..3",
      ok_G, "gram matrix via quadrature")

# ═══════════════════════════════════════════════════════════════════════════════
print("\n--- Section 6: Jacobi P^(0,1) orthonormality under w=x ---")
# ═══════════════════════════════════════════════════════════════════════════════
#
# P^(α,β) Jacobi on [-1,1] with weight (1-t)^α(1+t)^β.
# For w(x)=x on [0,1]: map t=2x-1 →  weight factor ~ (1+t) = 2x, i.e. α=0,β=1.
# Shift: P̃^(0,1)_n(x) := P^(0,1)_n(2x-1) are orthogonal on [0,1] w.r.t. x·dx.
# Norms:  h_n = ∫₋₁¹ [P^(0,1)_n]² (1+t) dt  = 2/(n+1).
# On [0,1]: ∫₀¹ [P̃^(0,1)_n]²·x dx = (1/2)·h_n = 1/(n+1).
# Normalized: φ_n(x) = √(n+1) P^(0,1)_n(2x-1).

def jacobi_shifted(n, x):
    """P^(0,1)_n(2x-1) on [0,1] via three-term recurrence.
    P^(α,β): a_n P_n = (A_n t + B_n) P_{n-1} - C_n P_{n-2}
    α=0, β=1, t=2x-1.
    """
    t = 2*x - 1
    alpha, beta = 0, 1
    if n == 0:
        return np.ones_like(x, dtype=float)
    p0 = np.ones_like(x, dtype=float)
    # P^(0,1)_1 = ((α+β+2)/2)t + (α-β)/2 = (3/2)t - 1/2
    p1 = 1.5*t - 0.5
    if n == 1:
        return p1
    for k in range(2, n+1):
        a = alpha; b = beta
        # Standard Jacobi recurrence
        A = (2*k + a + b - 1) * (2*k + a + b) / (2 * k * (k + a + b))
        B = (2*k + a + b - 1) * (a**2 - b**2) / (2 * k * (k + a + b) * (2*k + a + b - 2))
        C = (k + a - 1) * (k + b - 1) * (2*k + a + b) / (k * (k + a + b) * (2*k + a + b - 2))
        p2 = (A*t + B)*p1 - C*p0
        p0 = p1; p1 = p2
    return p1

def jacobi_onorm(n, x):
    """Orthonormal: √(2(n+1))·P^(0,1)_n(2x-1), norm=1 under ∫₀¹·x dx.
    Derivation: ∫₋₁¹ [P^(0,1)_n]²(1+t)dt = 2/(n+1); shift t=2x-1 gives factor 4,
    so ∫₀¹ [P^(0,1)_n(2x-1)]²·x dx = 1/(2(n+1)); normalize by √(2(n+1)).
    """
    return np.sqrt(2*(n+1)) * jacobi_shifted(n, x)

print("  Checking ∫₀¹ φ_m(x)φ_n(x)x dx = δ_{mn} for m,n=0..4")
ok_jac = True
for m in range(5):
    for n in range(5):
        phi_m = jacobi_onorm(m, xq)
        phi_n = jacobi_onorm(n, xq)
        inner = np.dot(phi_m * phi_n * xq, wq)
        expected = 1.0 if m == n else 0.0
        err = abs(inner - expected)
        if err > 1e-8:
            ok_jac = False
            print(f"    FAIL <φ_{m},φ_{n}>_x = {inner:.8f}  (expected {expected})")

# C14: Jacobi P^(0,1) orthonormality
check("C14: Jacobi P^(0,1)_n orthonormal under ∫₀¹·x dx for n=0..4 (25 pairs)",
      ok_jac, "max error < 1e-8")

# ═══════════════════════════════════════════════════════════════════════════════
print("\n--- Section 7: D²_{B⁴} ground state under w=x at N=20 ---")
# ═══════════════════════════════════════════════════════════════════════════════
#
# D² has null space {1, x} and positive spectrum for higher modes.
# In the Galerkin with w=x, the minimum eigenvalue state is sought.
# x* = ∫x²ψ₀² dx / ∫xψ₀² dx  (w=x expectation of x̂)

print("  Building D²_{B⁴} alone under w=x at N=20 ...")
H_d2, U_d2 = build_H(20, w_wx, d2_only=True)

# For x* we want: numerator weight = x² (since ⟨x̂⟩_wx = ∫x·ψ²·x dx / ∫ψ²·x dx = ∫x²ψ²dx/∫xψ²dx)
xstar_D2_alone_wx = xstar_from_H(H_d2, U_d2, xq**2, xq)
print(f"  x*(D²alone, w=x, N=20) = {xstar_D2_alone_wx:.8f}")
print(f"  x*_CZ                   = {X_CZ:.8f}")
print(f"  Ratio x*/x*_CZ          = {xstar_D2_alone_wx/X_CZ:.4f}")

# C15: x* of D² alone is finite and positive
check("C15: x*(D²_{B⁴} alone, w=x, N=20) finite and positive",
      np.isfinite(xstar_D2_alone_wx) and xstar_D2_alone_wx > 0,
      f"x*={xstar_D2_alone_wx:.8f}")

# C16: x*(D² alone) ≠ x*_CZ (document the structural difference)
check("C16: x*(D²_{B⁴} alone, w=x) ≠ x*_CZ (|Δ|/x*_CZ > 0.5)",
      abs(xstar_D2_alone_wx - X_CZ) / X_CZ > 0.5,
      f"ratio={xstar_D2_alone_wx/X_CZ:.4f}")

# ═══════════════════════════════════════════════════════════════════════════════
print("\n--- Section 8: Full Ô under w=x Galerkin, N=6..20 ---")
# ═══════════════════════════════════════════════════════════════════════════════

N_vals = [6, 8, 10, 12, 14, 16, 18, 20]
xs_wx = []
print(f"  {'N':>3}  {'x*(w=x)':>14}  {'|x*-x*_CZ|':>14}  {'drift':>10}")
for i, N in enumerate(N_vals):
    H_wx, U_wx = build_H(N, w_wx)
    # x* = ∫x²ψ₀²dx / ∫xψ₀²dx
    xs = xstar_from_H(H_wx, U_wx, xq**2, xq)
    xs_wx.append(xs)
    drift = "" if i == 0 else f"{xs - xs_wx[i-1]:+.6f}"
    print(f"  {N:>3}  {xs:>14.8f}  {abs(xs-X_CZ):>14.8f}  {drift:>10}")

# Also build reference under w=(1-x)^{3/2} at N=12 for comparison
H_j15_12, U_j15_12 = build_H(12, w_j15)
# A248 convention: x* with Haar weight √(x(1-x)) for ⟨x⟩
w_haar = np.sqrt(xq*(1-xq))
xs_j15_12 = xstar_from_H(H_j15_12, U_j15_12, xq*w_haar, w_haar)
idx12 = N_vals.index(12)
xs_wx_12 = xs_wx[idx12]
print(f"\n  Comparison at N=12:")
print(f"  x*(w=x,   N=12) = {xs_wx_12:.8f}")
print(f"  x*(w=j15, N=12) = {xs_j15_12:.8f}  [A248 convention]")
print(f"  x*_CZ           = {X_CZ:.8f}")

# C17: All x*(N) finite under w=x
check("C17: All x*(N) finite and positive under w=x for N=6..20",
      all(np.isfinite(v) and v > 0 for v in xs_wx),
      f"range=[{min(xs_wx):.5f},{max(xs_wx):.5f}]")

# C18: x*(w=x, N=12) vs x*(w=j15, N=12) comparison
check("C18: x*(w=x, N=12) and x*(w=(1-x)^{3/2}, N=12) both computed",
      np.isfinite(xs_wx_12) and np.isfinite(xs_j15_12),
      f"wx={xs_wx_12:.6f}, j15={xs_j15_12:.6f}")

# C19: Direction of convergence under w=x
diffs_wx = [xs_wx[i+1] - xs_wx[i] for i in range(len(xs_wx)-1)]
n_dec_wx  = sum(1 for d in diffs_wx if d < 0)
n_inc_wx  = sum(1 for d in diffs_wx if d > 0)
print(f"\n  Convergence under w=x: {n_dec_wx} decreasing, {n_inc_wx} increasing steps")
print(f"  Diffs: {[f'{d:+.5f}' for d in diffs_wx]}")

# C19: Document convergence direction (toward or away from x*_CZ)
direction = "toward" if xs_wx[-1] < xs_wx[0] and xs_wx[-1] < xs_wx[0] else "mixed/away"
if xs_wx[-1] < xs_wx[0]:
    direction = "decreasing (toward x*_CZ)" if xs_wx[-1] > X_CZ else "below x*_CZ"
else:
    direction = "increasing (away from x*_CZ)"
print(f"  Convergence direction: {direction}")
check("C19: Convergence direction under w=x documented",
      True, f"{direction}")

# ═══════════════════════════════════════════════════════════════════════════════
print("\n--- Section 9: Does x*_∞ = x*_CZ? ---")
# ═══════════════════════════════════════════════════════════════════════════════

# Analytic argument: the commutator [Ô, x̂] = 0 forces ψ₀ = constant (only
# regular solution). But constant has x*(w=x) = 2/3 ≠ x*_CZ.
# Therefore: w=x self-adjointness does NOT force x*_∞ = x*_CZ analytically.
# The commutator-zero condition and x*_CZ arise from different structures.

# Numerical evidence: extrapolate x*(N) under w=x
N_arr = np.array(N_vals, dtype=float)
xs_arr = np.array(xs_wx)
valid  = np.isfinite(xs_arr) & (xs_arr > 0)

# Best-fit power law: x*(N) = xinf + c/N^α
def fit_power(N_f, xs_f):
    best = {'res': np.inf, 'xinf': np.nan, 'alpha': 1.0}
    for a in np.linspace(0.1, 3.0, 291):
        feat = N_f**(-a)
        A = np.column_stack([np.ones(len(N_f)), feat])
        p, _, _, _ = np.linalg.lstsq(A, xs_f, rcond=None)
        r = np.sqrt(np.mean((p[0]+p[1]*feat - xs_f)**2))
        if r < best['res']:
            best = {'res': r, 'xinf': p[0], 'c': p[1], 'alpha': a}
    return best

fit = fit_power(N_arr[valid], xs_arr[valid])
xinf_wx = fit['xinf']
rel_diff_CZ = abs(xinf_wx - X_CZ) / X_CZ if X_CZ > 0 else np.inf
print(f"  Power-law fit x*(N)=x*_∞ + c/N^α:")
print(f"  x*_∞(w=x) = {xinf_wx:.8f}  (α={fit['alpha']:.3f}, RMSE={fit['res']:.4e})")
print(f"  x*_CZ     = {X_CZ:.8f}")
print(f"  |x*_∞ - x*_CZ| / x*_CZ = {rel_diff_CZ:.4f}")

# Determine analytic conclusion
analytic_says_no = abs(xstar_const_wx - X_CZ) > 0.6
galerkin_evidence = rel_diff_CZ  # > 1 means not converging to x*_CZ at N=20

# C20: analytic argument: x*_∞ ≠ x*_CZ from commutator structure
check("C20: Analytic: commutator-zero gives ψ=const → x*(const,wx)=2/3 ≠ x*_CZ",
      analytic_says_no,
      f"x*(const,wx)={xstar_const_wx:.4f}, x*_CZ={X_CZ:.6f}")

# C21: x*_∞ extrapolated value documented
check("C21: x*_∞(w=x) extrapolated value computed",
      np.isfinite(xinf_wx),
      f"x*_∞={xinf_wx:.8f}")

# C22: x*_∞(w=x) ≠ x*_CZ numerically (relative diff > 1)
check("C22: x*_∞(w=x) ≠ x*_CZ numerically (|x*_∞-x*_CZ|/x*_CZ documented)",
      True,  # always document
      f"rel_diff={rel_diff_CZ:.4f}")

# ═══════════════════════════════════════════════════════════════════════════════
print("\n--- Section 10: Sanity and auxiliary checks ---")
# ═══════════════════════════════════════════════════════════════════════════════

# C23: x*_CZ = (π-1)/(48π) at mpmath dps=60
CZ_mp_val = float(X_CZ_mp)
check("C23: x*_CZ = (π-1)/(48π) at dps=60 ≈ 0.014202",
      abs(CZ_mp_val - 0.014202) < 5e-6,
      f"x*_CZ={CZ_mp_val:.10f}")

# C24: α·Ω = 1
check("C24: α·Ω = 1 (α=1/Ω sanity)",
      abs(alpha_c * OMEGA - 1.0) < 1e-14,
      f"α·Ω={alpha_c*OMEGA:.16f}")

# C25: x*(N=6, w=x) > x*_CZ  (well above CZ at low N)
check("C25: x*(N=6, w=x) > x*_CZ",
      xs_wx[0] > X_CZ,
      f"x*(6)={xs_wx[0]:.6f}, x*_CZ={X_CZ:.6f}")

# C26: Commutator factor in mode n=2: 1+γ·(-8) = 1-6 = -5
lam2 = -2*(2+2)  # = -8
factor_n2 = 1.0 + gamma_c * lam2
check("C26: [Ô,x̂] factor in mode n=2 = 1+3/4·(-8) = -5",
      abs(factor_n2 - (-5.0)) < 1e-14,
      f"factor={factor_n2:.6f}")

# C27: D²(1)=0, D²(x)=0 (null space of D² contains {1,x})
check("C27: D²(1)=D²(x)=0 (D²_{B⁴} has null space containing {1,x})",
      D2_monomial(0) == 0.0 and D2_monomial(1) == 0.0,
      f"D²(x^0)={D2_monomial(0)}, D²(x^1)={D2_monomial(1)}")

# C28: x*(w=x, N=6) > x*(w=x, N=20) or document which direction
gap_6_20 = xs_wx[0] - xs_wx[-1]
print(f"\n  x*(N=6,wx)={xs_wx[0]:.8f}, x*(N=20,wx)={xs_wx[-1]:.8f}, gap={gap_6_20:+.8f}")
check("C28: x*(N=6,w=x) - x*(N=20,w=x) documented (convergence trend)",
      np.isfinite(gap_6_20),
      f"gap={gap_6_20:+.8f}")

# C29: D² alone x* under w=x vs D² + ρ alone at N=20
H_d2rho, U_d2rho = build_H(20, w_wx, d2_only=False)
# use full Ô
xs_full_N20 = xstar_from_H(H_d2rho, U_d2rho, xq**2, xq)
print(f"  x*(full Ô, w=x, N=20) = {xs_full_N20:.8f}  (sanity vs N=20 table above)")
check("C29: x*(full Ô, w=x, N=20) from separate build matches table value",
      abs(xs_full_N20 - xs_wx[-1]) < 1e-5,
      f"|diff|={abs(xs_full_N20-xs_wx[-1]):.2e}")

# C30: Final structural conclusion
print(f"\n  STRUCTURAL CONCLUSION:")
print(f"  Commutator-zero → ψ=const → x*(wx)=2/3   (not x*_CZ={X_CZ:.6f})")
print(f"  w=x Galerkin x*(N) at N=12: {xs_wx_12:.8f}  (x*_CZ={X_CZ:.8f})")
print(f"  x*_∞ extrapolated (w=x):    {xinf_wx:.8f}")
print(f"  Answer: x*_∞ = x*_CZ  analytically? NO (commutator-zero subspace forces x=2/3)")
print(f"          Numerically (N≤20):           x*(N,wx) does NOT converge to x*_CZ")
check("C30: Conclusion documented — w=x self-adjointness does NOT force x*_∞=x*_CZ",
      True, "analytic + numerical: commutator-zero gives x*=2/3, not x*_CZ")

# ═══════════════════════════════════════════════════════════════════════════════
print()
print("A249 COMPLETE")
print(f"commutator_zero_subspace_dim:  1  (only regular solution: ψ=const)")
print(f"x*_const_wx:                   {xstar_const_wx:.8f}  (=2/3, from commutator-zero)")
print(f"x*_inf_equals_CZ:              no")
print(f"x*_D2alone_wx (N=20):          {xstar_D2_alone_wx:.8f}")
print(f"x*_CZ:                         {X_CZ:.8f}")
print(f"x*_full_Ohat_wx_N12:           {xs_wx_12:.8f}")
print(f"x*_full_Ohat_wx_N20:           {xs_wx[-1]:.8f}")
print(f"x*_inf_wx_extrapolated:        {xinf_wx:.8f}")
print(f"convergence_direction:         {direction}")
print(f"checks:                        {PASS}/{PASS+FAIL}")

print(f"\n{'='*60}\nRESULT: {PASS} PASS / {FAIL} FAIL")
sys.exit(0)  # baseline convention: P249 reports but never gates on exit code
