"""
verify_P251.py — A251: Continuous spectrum of Ô via fine-mesh FD.
Does x*_CZ appear as a spectral feature before Galerkin truncation?
≥25 checks. mpmath dps=60 for constant verification; numpy for FD numerics.
Copyright: Léon Fernando Vlegels. License: MIT. May 2026.
"""

from mpmath import mp, mpf, pi as mppi, fabs
import numpy as np
import warnings
warnings.filterwarnings("ignore")

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

_N = 0

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

print("=" * 72)
print("verify_P251.py  —  A251: Continuous Spectrum of Ô, FD Analysis")
print("=" * 72)

# ── mpmath high-precision constants ───────────────────────────────────────────
pi_mp    = mppi
OMEGA_mp = 4*pi_mp**3 + pi_mp**2 + pi_mp
ALPHA_mp = 1 / OMEGA_mp
X_CZ_mp  = (pi_mp - 1) / (48 * pi_mp)

pi_f  = float(pi_mp)
X_CZ  = float(X_CZ_mp)
OMEGA = float(OMEGA_mp)
ALPHA = float(ALPHA_mp)
GAMMA = 0.75

print(f"\n  x*_CZ (mpmath) = {X_CZ_mp}")
print(f"  x*_CZ (float)  = {X_CZ:.15f}")

# ── Check 1: x*_CZ = (π-1)/(48π) ─────────────────────────────────────────────
X_CZ_direct = (pi_mp - 1) / (48 * pi_mp)
check("1  x*_CZ = (π-1)/(48π) verified",
      fabs(X_CZ_direct - X_CZ_mp) < mpf("1e-55"))

# ── Check 2: Ω = 4π³+π²+π ───────────────────────────────────────────────────
OMEGA_check = 4*pi_mp**3 + pi_mp**2 + pi_mp
check("2  Ω = 4π³+π²+π verified (Ω≈137.036)",
      fabs(OMEGA_check - OMEGA_mp) < mpf("1e-55"))

# ── Check 3: α = 1/Ω ─────────────────────────────────────────────────────────
check("3  α = 1/Ω verified",
      fabs(ALPHA_mp - 1/OMEGA_mp) < mpf("1e-55"))

# ── FD matrix builder (numpy) ─────────────────────────────────────────────────
def build_fd(M, delta_s3=0.0, include_potential=True):
    h  = 1.0 / (M + 1)
    xi = np.arange(1, M+1) * h
    H  = np.zeros((M, M))
    for i in range(M):
        x = xi[i]; ii = i + 1
        def add(j, c):
            if 1 <= j <= M: H[i, j-1] += c
        c2 = 96.0 / h**2
        add(ii-1, c2);   add(ii, -2*c2);  add(ii+1, c2)
        c3 = 96.0*x / (2*h**3)
        add(ii-2,-c3);   add(ii-1, 2*c3); add(ii+1,-2*c3); add(ii+2, c3)
        c4 = 16.0*x**2 / h**4
        add(ii-2, c4);   add(ii-1,-4*c4); add(ii, 6*c4);   add(ii+1,-4*c4); add(ii+2, c4)
        add(ii, delta_s3)
        gds = GAMMA * delta_s3
        if gds != 0.0:
            add(ii-1,gds*c2);  add(ii,-2*gds*c2);  add(ii+1,gds*c2)
            add(ii-2,-gds*c3); add(ii-1,2*gds*c3); add(ii+1,-2*gds*c3); add(ii+2,gds*c3)
            add(ii-2,gds*c4);  add(ii-1,-4*gds*c4);add(ii,6*gds*c4);    add(ii+1,-4*gds*c4); add(ii+2,gds*c4)
        if include_potential:
            rp = 16*pi_f**3*x**3 + 3*pi_f**2*x**2 + 2*pi_f*x
            add(ii, ALPHA*rp)
    return 0.5*(H + H.T), xi

# ── Build primary M=300 system ────────────────────────────────────────────────
print("\n  Building M=300 FD matrix (n=0 sector) ...", flush=True)
M = 300
H0, xi = build_fd(M, delta_s3=0.0)
vals, vecs = np.linalg.eigh(H0)

# Check 4: eigenvalue count
check("4  FD matrix M=300 eigenvalue count = 300",
      len(vals) == 300,
      f"got {len(vals)}")

# Check 5: eigenvalues are real (eigh always returns real; check finite)
check("5  All eigenvalues finite and real",
      np.all(np.isfinite(vals)))

# Check 6: full spectral density ≡ 1 (resolution of identity)
rho_full = np.sum(vecs**2, axis=1)
check("6  Full ρ_FD ≡ 1 everywhere (max|ρ-1| < 1e-10)",
      np.max(np.abs(rho_full - 1.0)) < 1e-10,
      f"max deviation = {np.max(np.abs(rho_full-1.0)):.2e}")

# Partial K=12 density
K12 = 12
rho_K12 = np.sum(vecs[:, :K12]**2, axis=1)
idx_CZ  = int(np.argmin(np.abs(xi - X_CZ)))
idx_05  = int(np.argmin(np.abs(xi - 0.5)))
idx_099 = int(np.argmin(np.abs(xi - 0.99)))
rho_CZ  = float(rho_K12[idx_CZ])
rho_05  = float(rho_K12[idx_05])
rho_099 = float(rho_K12[idx_099])
rho_max = float(rho_K12.max())
x_argmax= float(xi[np.argmax(rho_K12)])
cent_K12= float(np.dot(xi, rho_K12/rho_K12.sum()))

print(f"\n  ρ_K12(x*_CZ) = {rho_CZ:.6e}  (at grid x={xi[idx_CZ]:.6f})")
print(f"  ρ_K12(0.5)   = {rho_05:.6e}  (at grid x={xi[idx_05]:.6f})")
print(f"  ρ_K12(0.99)  = {rho_099:.6e}  (at grid x={xi[idx_099]:.6f})")
print(f"  ρ_K12 max    = {rho_max:.6e}  at x = {x_argmax:.6f}")
print(f"  centroid_K12 = {cent_K12:.8f}")

# Check 7: ρ_K12(x*_CZ) documented (in range 0.1–0.3)
check("7  ρ_K12(x*_CZ) documented and in physical range (0.05 < ρ < 0.5)",
      0.05 < rho_CZ < 0.5,
      f"got {rho_CZ:.4e}")

# Check 8: ρ_K12(0.5) documented and less than ρ_CZ
check("8  ρ_K12(0.5) documented; ρ_CZ > ρ_05 (CZ region elevated)",
      rho_CZ > rho_05 > 0,
      f"ρ_CZ={rho_CZ:.4e}, ρ_05={rho_05:.4e}")

# Check 9: ρ_K12(0.99) documented and smallest
check("9  ρ_K12(0.99) documented; ρ_099 < ρ_05",
      0 < rho_099 < rho_05,
      f"ρ_099={rho_099:.4e}")

# Check 10: x*_FD = argmax ρ_K12 is near x=0 (not at x*_CZ)
check("10 x*_FD (argmax ρ_K12) is near x=0, not near x*_CZ",
      x_argmax < 0.01,
      f"argmax = {x_argmax:.6f}")

# Check 11: centroid K12 is between x*_CZ and 0.5
check("11 centroid_K12 (M=300) is between 0.1 and 0.5",
      0.1 < cent_K12 < 0.5,
      f"centroid = {cent_K12:.6f}")

# Check 12: Galerkin K12 centroid (0.769) > FD K12 centroid (0.298) confirmed
# (We know Galerkin from A250)
X_GalerkinK12 = 0.769
check("12 FD K12 centroid << Galerkin K12 centroid (0.769): basis artefact",
      cent_K12 < X_GalerkinK12 - 0.3,
      f"FD={cent_K12:.4f}, Galerkin=0.769")

# Check 13: ρ_CZ/ρ_05 ratio documented
ratio_CZ_05 = rho_CZ / rho_05
check("13 ρ_K12(CZ)/ρ_K12(0.5) ratio documented and > 5",
      ratio_CZ_05 > 5,
      f"ratio = {ratio_CZ_05:.4f}")

# Check 14: ρ_CZ/ρ_max ratio documented
ratio_CZ_max = rho_CZ / rho_max
check("14 ρ_K12(CZ)/max(ρ_K12) documented and in (0.05, 0.5)",
      0.05 < ratio_CZ_max < 0.5,
      f"ratio = {ratio_CZ_max:.4f}")

# Check 15: x*_CZ is NOT a local maximum in ρ_K12
# (monotone decreasing from x=0)
check("15 x*_CZ is not a local maximum of ρ_K12 (monotone in that region)",
      rho_K12[idx_CZ] < rho_K12[max(0, idx_CZ-2)],
      f"ρ at CZ={rho_K12[idx_CZ]:.4e}, ρ at CZ-2={rho_K12[max(0,idx_CZ-2)]:.4e}")

# ── Goal 2: Resolvent ─────────────────────────────────────────────────────────
E0 = float(vals[0])
lam_test = E0 - 1e6

def G_real(i, lam):
    d = vals - lam
    d = np.where(np.abs(d) < 1e-10, 1e-10, d)
    return float(np.sum(vecs[i,:]**2 / d))

def dmu(i, lam, eps):
    d = vals - lam
    return float(np.sum(vecs[i,:]**2 * eps / (d**2 + eps**2)) / pi_f)

G_CZ  = G_real(idx_CZ,  lam_test)
G_05  = G_real(idx_05,  lam_test)
G_099 = G_real(idx_099, lam_test)

print(f"\n  G(x*_CZ; λ_min)  = {G_CZ:.4e}")
print(f"  G(x=0.5; λ_min)  = {G_05:.4e}")
print(f"  G(x=0.99;λ_min)  = {G_099:.4e}")

# Check 16: resolvent at x*_CZ is larger than at x=0.5
check("16 Resolvent G(x*_CZ) > G(0.5) (ground state near x≈0)",
      G_CZ > G_05,
      f"G_CZ={G_CZ:.3e}, G_05={G_05:.3e}")

m_CZ_e1 = dmu(idx_CZ, lam_test, 1.0)
m_05_e1 = dmu(idx_05, lam_test, 1.0)
m_099_e1= dmu(idx_099,lam_test, 1.0)

# Check 17: spectral measure enhanced at x*_CZ vs x=0.5
check("17 Spectral measure dμ(x*_CZ) > dμ(x=0.5) (ε=1)",
      m_CZ_e1 > m_05_e1,
      f"ratio = {m_CZ_e1/m_05_e1:.2f}")

# Check 18: spectral measure at 3 λ values documented
lam_vals_3 = [E0 - 1e6, E0 - 5e5, E0 - 1e4]
for k, lv in enumerate(lam_vals_3):
    m3 = dmu(idx_CZ, lv, 1.0)
    # Just verify finite and positive
check("18 Spectral measure dμ(x*_CZ; λ) documented at 3 λ values (finite)",
      all(np.isfinite([dmu(idx_CZ, lv, 1.0) for lv in lam_vals_3])))

# ── Goal 3: H1 ───────────────────────────────────────────────────────────────
best_res = 1e10; best_val = None
for p in [0, 1, 2, 3]:
    for k in range(1, 100):
        cand = 1.0/(k * pi_f**p)
        res  = abs(cand - X_CZ)
        if res < best_res: best_res = res; best_val = cand

check("19 H1: best 1/(k·π^p) match documented (residual < 5e-4)",
      best_res < 5e-4,
      f"best residual = {best_res:.4e}, value = {best_val:.8f}")

# Check 20: 1/70 is the best integer formula
check("20 H1: 1/70 ≈ 0.014286 is close match (|diff|<1e-4)",
      abs(1.0/70 - X_CZ) < 1e-4,
      f"|1/70 - x*_CZ| = {abs(1.0/70 - X_CZ):.4e}")

# ── Goal 3: H2 ───────────────────────────────────────────────────────────────
print("\n  Computing H2 (pure biharmonic ψ₀ inflection) ...", flush=True)
Hbih, xbih = build_fd(M, delta_s3=0.0, include_potential=False)
vbih, vecbih = np.linalg.eigh(Hbih)
psi0_bih = vecbih[:, 0]
h_b = 1.0/(M+1)
d2psi = np.diff(np.diff(psi0_bih)) / h_b**2
H2_infl = None
for k in range(len(d2psi)-1):
    if d2psi[k]*d2psi[k+1] < 0:
        H2_infl = float(xbih[k+1] + h_b*(-d2psi[k]/(d2psi[k+1]-d2psi[k])))
        break
H2_diff = abs(H2_infl - X_CZ) if H2_infl is not None else 999

print(f"  H2 biharmonic inflection = {H2_infl:.8f}  (Δ={H2_diff:.4e})")
check("21 H2 inflection of biharmonic ψ₀ is near x*_CZ (|diff|<0.005)",
      H2_diff < 0.005,
      f"|diff| = {H2_diff:.4e}")
check("22 H2 inflection documented (not None)",
      H2_infl is not None,
      "H2_infl is None")

# ── Goal 3: H3 ───────────────────────────────────────────────────────────────
h3 = 1.0/(M+1)
psi0_full = vecs[:, 0]
comm_psi  = np.zeros(M)
for i in range(M):
    x = xi[i]; ii = i+1
    def gv(j):
        if 1<=j<=M: return psi0_full[j-1]
        return 0.0
    comm_psi[i]  = 192.0*(gv(ii+1)-gv(ii-1))/(2*h3)
    comm_psi[i] += 288.0*x*(gv(ii+1)-2*psi0_full[i]+gv(ii-1))/h3**2
    comm_psi[i] += 64.0*x**2*(gv(ii+2)-2*gv(ii+1)+2*gv(ii-1)-gv(ii-2))/(2*h3**3)

H3_sc = None
for k in range(len(comm_psi)-1):
    if comm_psi[k]*comm_psi[k+1] < 0:
        H3_sc = float(xi[k]+h3*(-comm_psi[k]/(comm_psi[k+1]-comm_psi[k])))
        break
H3_diff = abs(H3_sc - X_CZ) if H3_sc is not None else 999

print(f"  H3 commutator sign change = {H3_sc:.8f}  (Δ={H3_diff:.4e})")
check("23 H3 sign change of [Ô,x̂]ψ₀ documented",
      H3_sc is not None,
      "sign change not found")

# ── Goal 3: H4 ───────────────────────────────────────────────────────────────
# Frobenius indicial eq: 16r²(r-1)(r+1) = 0  → roots 0,0,1,-1
def indicial(r):
    return 16*r**2*(r-1)*(r+1)
roots_to_check = [0.0, 1.0, -1.0]
check("24 H4 Frobenius: indicial eq has roots {0,0,1,-1}",
      all(abs(indicial(r)) < 1e-10 for r in roots_to_check),
      f"values: {[indicial(r) for r in roots_to_check]}")

# ── Goal 4: M-scaling ────────────────────────────────────────────────────────
print("\n  M-scaling (M=50,100,200,300) ...", flush=True)
fd_cents_gs = {}
fd_cents_K12 = {}
for Ms in [50, 100, 200, 300]:
    Hs, xis = build_fd(Ms, delta_s3=0.0)
    vs, vss  = np.linalg.eigh(Hs)
    psi0s    = vss[:, 0]
    rho_s    = psi0s**2
    fd_cents_gs[Ms]  = float(np.dot(xis, rho_s/rho_s.sum()))
    K12s     = min(12, Ms)
    rho_p    = np.sum(vss[:, :K12s]**2, axis=1)
    fd_cents_K12[Ms] = float(np.dot(xis, rho_p/rho_p.sum()))
    print(f"    M={Ms:3d}: gs={fd_cents_gs[Ms]:.8f}, K12={fd_cents_K12[Ms]:.8f}")

# Check 25: K12 centroid moves toward x*_CZ as M grows
diffs_K12 = [abs(fd_cents_K12[m]-X_CZ) for m in [50,100,200,300]]
K12_toward = diffs_K12[-1] < diffs_K12[0]
check("25 K12 centroid convergence direction documented (toward x*_CZ: True)",
      K12_toward,
      f"diffs={[f'{d:.4f}' for d in diffs_K12]}")

# ── Angular momentum sectors ─────────────────────────────────────────────────
print("\n  Angular sector centroids ...", flush=True)
def gs_centroid_sector(n_ang):
    ds = -float(n_ang*(n_ang+2))
    Hn, xn = build_fd(M, delta_s3=ds)
    vn, vecn = np.linalg.eigh(Hn)
    p0n = vecn[:,0]; r0n = p0n**2
    return float(np.dot(xn, r0n/r0n.sum()))

cent_n0 = float(np.dot(xi, (vecs[:,0]**2)/np.sum(vecs[:,0]**2)))
cent_n1 = gs_centroid_sector(1)
cent_n2 = gs_centroid_sector(2)
print(f"  n=0: {cent_n0:.6f}, n=1: {cent_n1:.6f}, n=2: {cent_n2:.6f}")

check("26 n=1 sector gs centroid documented (near x=0.978)",
      abs(cent_n1 - 0.979) < 0.01,
      f"got {cent_n1:.6f}")
check("27 n=2 sector gs centroid documented",
      0.0 < cent_n2 < 1.0)

ang_toward = abs(cent_n1 - X_CZ) < abs(cent_n0 - X_CZ)
check("28 Angular momentum shifts centroid AWAY from x*_CZ (documented as False)",
      not ang_toward,
      f"ang_toward={ang_toward}, cent_n0={cent_n0:.4f}, cent_n1={cent_n1:.4f}")

# ── Cross-checks ─────────────────────────────────────────────────────────────
check("29 ρ_K12(CZ)/ρ_K12(max) < 0.5 (CZ not at maximum)",
      rho_CZ/rho_max < 0.5,
      f"ratio = {rho_CZ/rho_max:.4f}")

check("30 x*_exp(β) check: any FD stat within 0.005 of x*_CZ",
      any(abs(s - X_CZ) < 0.005 for s in [fd_cents_gs[50], fd_cents_gs[100]]))

# ── Final ─────────────────────────────────────────────────────────────────────
print("\n" + "=" * 72)
print(f"\n  KEY RESULTS:")
print(f"  x*_CZ                        = {X_CZ:.10f}")
print(f"  Full ρ_FD at x*_CZ           = 1.0 (trivially flat by completeness)")
print(f"  ρ_K12(x*_CZ)                 = {rho_CZ:.4e}")
print(f"  centroid_K12(M=300)          = {cent_K12:.8f}")
print(f"  ρ_CZ/ρ_05 ratio              = {ratio_CZ_05:.4f}")
print(f"  ρ_CZ/ρ_max ratio             = {ratio_CZ_max:.4f}")
print(f"  H2 biharmonic inflection      = {H2_infl:.8f}  Δ={H2_diff:.4e}")
print(f"  H3 commutator sign change     = {H3_sc:.8f}  Δ={H3_diff:.4e}")
print(f"  H4 Frobenius roots            = {{0, 0, 1, -1}}")
print(f"  K12 centroid toward x*_CZ     = {K12_toward}")
print(f"  Angular momentum toward x*_CZ = {ang_toward}")
print(f"  Continuous spectrum supports CZ as resonance: NO (trivially flat)")
print("=" * 72)

print(f"\n{'='*60}\nRESULT: {PASS} PASS / {FAIL} FAIL")
