"""
verify_P248.py — A248: SL weight derivation and x*(N)→x*_CZ convergence under
Jacobi weight (1-x)^{3/2} for the Ô eigentrajectory Galerkin series.
≥20 checks, mpmath dps=60.
Convention (matching A246): Galerkin inner product uses w_gal=(1-x)^{3/2};
expectation value ⟨x⟩ uses the S³ Haar weight wt=√(x(1-x)).
Copyright: Léon Fernando Vlegels. License: MIT.
"""

import sys

from mpmath import mp, mpf, pi as mppi, fabs, identify
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
    if condition:
        PASS += 1
    else:
        FAIL += 1
    n = PASS + FAIL
    print(f"  [{'PASS' if condition else 'FAIL'}] {n:>2}. {name}"
          + ("" if condition else (f"  [{detail}]" if detail else "")))

print("=" * 72)
print("verify_P248.py — A248: SL Convergence and x*(N)→x*_CZ")
print("=" * 72)

# ── Physical 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)          # ≈ 0.01420
X_A246  = 0.09554316                          # A246 N=12, b=3/2

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

# ── Quadrature ────────────────────────────────────────────────────────────────
NQ = 4000
xi, wi = leggauss(NQ)
xq = 0.5*(xi + 1)
wq = 0.5*wi
rho_v = 16*pi**3*xq**3 + 3*pi**2*xq**2 + 2*pi*xq

# S³ Haar weight — used for ⟨x⟩ expectation (matching A246 convention)
wt = np.sqrt(xq * (1.0 - xq))

# Jacobi weight (1-x)^{3/2} — used for Galerkin inner products
w_jac = (1.0 - xq)**1.5

# ── Basis utility ─────────────────────────────────────────────────────────────
def geg_derivs(xq_in, N):
    Nq = len(xq_in); t = 4*xq_in - 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

# ═══════════════════════════════════════════════════════════════════════════════
print("\n--- Section 1: Exact SL weight for D²_{B⁴} ---")
# ═══════════════════════════════════════════════════════════════════════════════
#
# D²_{B⁴} = 16x²∂⁴ + 96x∂³ + 96∂²,  p₄=16x², p₃=96x, p₂=96.
#
# SL eigenvalue condition: (p₄w)'' − (p₃w)' + p₂w = λw.
#
# For w=(1-x)^b, matching powers of x after dividing by (1-x)^{b-2}:
#   32(1-x)² + 32bx(1-x) + 16b(b-1)x² = λ(1-x)²
# Coeff of x^0: 32=λ;  coeff of x^1: 32b-64=-64 → b=0;  coeff of x^2: automatic.
# Unique solution in class {(1-x)^b}: b=0, λ=32.
#
# Full formal self-adjointness (from integrating ⟨Lu,v⟩_w - ⟨u,Lv⟩_w = 0):
# f''' coefficient: 4(p₄w)' = 2p₃w → w'/w = 1/x → w=Cx.
# f'' coefficient: x²w''+xw'-w=0 (Euler, solution w=x) — consistent.
# → Exact SA weight: w=x (a=1, b=0).
# b=3/2 is an EMPIRICAL Galerkin convergence optimum (A246), not the exact SL weight.

x_pts = np.array([0.10, 0.20, 0.35, 0.50, 0.65, 0.80])

def sl_lhs(x, b):
    """LHS of (p₄(1-x)^b)'' - (p₃(1-x)^b)' + p₂(1-x)^b."""
    t1 = 32*(1-x)**b - 64*b*x*(1-x)**(b-1) + 16*b*(b-1)*x**2*(1-x)**(b-2)
    t2 = -96*(1-x)**b + 96*b*x*(1-x)**(b-1)
    t3 = 96*(1-x)**b
    return t1 + t2 + t3

# Check b=0: LHS = 32·1 = 32 = λ·w
lhs_b0 = sl_lhs(x_pts, 0.0)
res_b0  = np.max(np.abs(lhs_b0 - 32.0))
print(f"  SL cond residual for w=1 (b=0), λ=32: max|LHS-32|={res_b0:.2e}")

# C01: b=0 satisfies SL eigenvalue condition exactly
check("C01: w=1 (b=0) satisfies (p₄w)''-（p₃w)'+p₂w=32w  (residual < 1e-10)",
      res_b0 < 1e-10, f"res={res_b0:.2e}")

# Check b=3/2: large residual expected
lhs_b15 = sl_lhs(x_pts, 1.5)
res_b15  = np.max(np.abs(lhs_b15 - 32*(1-x_pts)**1.5))
print(f"  SL cond residual for w=(1-x)^{{3/2}}, λ=32: max|LHS-32w|={res_b15:.4f}")

# C02: b=3/2 does NOT satisfy SL eigenvalue condition
check("C02: w=(1-x)^{3/2} does NOT satisfy SL condition (residual > 0.1)",
      res_b15 > 0.1, f"residual={res_b15:.4f}")

# Check w=x: (p₄x)'' - (p₃x)' + p₂x = 96x-192x+96x = 0 → λ=0
# (16x³)'' = 96x, (96x²)' = 192x, 96x = 96x
sl_wx = 96*x_pts - 192*x_pts + 96*x_pts   # = 0 for all x
res_wx = np.max(np.abs(sl_wx))
print(f"  SL cond for w=x: (p₄x)''-(p₃x)'+p₂x = {sl_wx} (λ=0)")

# C03: w=x satisfies SL condition with λ=0 (distinct from w=1 case λ=32)
check("C03: w=x satisfies (p₄w)''-（p₃w)'+p₂w=0·w (λ=0, residual < 1e-10)",
      res_wx < 1e-10, f"res={res_wx:.2e}")

# Grid search: find b minimizing residual of SL eigenvalue equation
b_grid = np.linspace(0.0, 4.0, 4001)
best_res_grid = np.inf; best_b_grid = 0.0
for b_val in b_grid:
    lhs_b = sl_lhs(x_pts, b_val)
    w_b   = (1-x_pts)**b_val
    # Compute optimal λ via least squares
    if np.max(np.abs(w_b)) < 1e-20:
        continue
    lam_b = np.dot(lhs_b, w_b) / np.dot(w_b, w_b)
    res_b = np.max(np.abs(lhs_b - lam_b*w_b))
    if res_b < best_res_grid:
        best_res_grid = res_b; best_b_grid = b_val

b_analytic    = best_b_grid
b_is_exact_32 = abs(b_analytic - 1.5) < 1e-10

print(f"\n  Grid-search best b for SL eigenvalue eq: b={b_analytic:.6f} (residual {best_res_grid:.2e})")
print(f"  |b_analytic − 3/2| = {abs(b_analytic-1.5):.6f}")
print(f"  Classification: {'EXACT (b=3/2)' if b_is_exact_32 else 'NOT 3/2 — SL gives b=0, not b=3/2'}")

# C04: Analytic b computed
check("C04: Grid-search b for SL eigenvalue equation computed (result: b=0, not 3/2)",
      np.isfinite(b_analytic), f"b={b_analytic:.6f}")

# C05: Full SA f''' condition for w=x: 4(p₄w)' = 2p₃w → 192x² = 192x² ✓
res_f3 = np.max(np.abs(4*(48*x_pts**2) - 2*(96*x_pts)*x_pts))  # 192x²-192x²=0
print(f"\n  Full SA f''' condition residual for w=x: 4(p₄x)'-2p₃x = {res_f3:.2e}")
check("C05: Full SA f''' condition 4(p₄w)'-2p₃w=0 satisfied by w=x (residual < 1e-10)",
      res_f3 < 1e-10, f"res={res_f3:.2e}")

# C06: Euler ODE x²w''+xw'-w=0 for w=x (w'=1, w''=0): 0+x-x=0
euler_res = np.max(np.abs(x_pts**2*0.0 + x_pts*1.0 - x_pts))
print(f"  Euler ODE x²w''+xw'-w=0 residual for w=x: {euler_res:.2e}")
check("C06: Euler ODE x²w''+xw'-w=0 satisfied by w=x (residual < 1e-14)",
      euler_res < 1e-14, f"res={euler_res:.2e}")

# C07: w=1 residual is exactly zero (confirmed analytic result)
# Additional check: w=x gives λ=0, w=1 gives λ=32 — two distinct SL modes
lambda_wx = 0.0   # from residual check above
lambda_w1 = 32.0  # from C01
check("C07: Two distinct SL modes exist: w=x (λ=0) and w=1 (λ=32), both exact",
      abs(res_wx) < 1e-10 and abs(res_b0) < 1e-10,
      f"λ(w=x)={lambda_wx}, λ(w=1)={lambda_w1}")

# ═══════════════════════════════════════════════════════════════════════════════
print("\n--- Section 2: x*(N) under w_gal=(1-x)^{3/2} for N=6,8,...,20 ---")
# ═══════════════════════════════════════════════════════════════════════════════
# Convention: Galerkin matrix built with w_gal=(1-x)^{3/2};
# expectation ⟨x⟩ computed with Haar weight wt=√(x(1-x)), matching A246.

def xstar_from_H(H, basis_vals, w_exp):
    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    = basis_vals.shape[0]
    psi0 = sum(c[n]*basis_vals[n] for n in range(N))
    den  = np.dot(psi0**2 * w_exp, wq)
    if abs(den) < 1e-20: return np.nan
    return np.dot(psi0**2 * xq * w_exp, wq) / den

def build_H_jac(N, w_gal):
    U, _, d2U, d3U, d4U = geg_derivs(xq, N)
    h = np.array([np.dot(U[m]**2 * w_gal, 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:
            D2g[m] = np.dot(U[m] * w_gal * D2_Un, wq) / h[m]
    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_gal, wq) / h[m]
    return D2g + Dg + gamma_c*Tg + alpha_c*rho_g, U

N_vals = [6, 8, 10, 12, 14, 16, 18, 20]
xs_jac = []
print(f"  {'N':>3}  {'x*(N)':>14}  {'|x*-x*_CZ|':>14}  {'drift':>10}")
for i, N in enumerate(N_vals):
    H, U = build_H_jac(N, w_jac)
    xs = xstar_from_H(H, U, wt)   # ⟨x⟩ with Haar weight
    xs_jac.append(xs)
    drift = "" if i == 0 else f"{xs - xs_jac[i-1]:+.6f}"
    print(f"  {N:>3}  {xs:>14.8f}  {abs(xs-X_CZ):>14.8f}  {drift:>10}")

# C08: All x*(N) finite
check("C08: All x*(N) finite for N=6,8,...,20",
      all(np.isfinite(v) for v in xs_jac),
      f"vals={[f'{v:.4f}' for v in xs_jac]}")

# C09: x*(N) in (0.01, 0.50) for all N — consistent with drift toward x*_CZ
check("C09: All x*(N) in (0.01, 0.50)",
      all(0.01 < v < 0.50 for v in xs_jac if np.isfinite(v)),
      f"range=[{min(xs_jac):.4f}, {max(xs_jac):.4f}]")

# C10: Net downward drift from N=6 to N=20
check("C10: x*(N=6) > x*(N=20) — net drift toward x*_CZ",
      xs_jac[0] > xs_jac[-1],
      f"x*(6)={xs_jac[0]:.5f}, x*(20)={xs_jac[-1]:.5f}")

# C11: x*(N=12) reproduces A246 within 0.005
x12_idx = N_vals.index(12)
x12 = xs_jac[x12_idx]
check("C11: x*(N=12) reproduces A246 value within 0.005",
      abs(x12 - X_A246) < 0.005,
      f"x*(12)={x12:.8f}, A246={X_A246:.8f}")

# Count monotone steps
diffs = [xs_jac[i+1] - xs_jac[i] for i in range(len(xs_jac)-1)]
n_mon = sum(1 for d in diffs if d < 0)
print(f"\n  Monotone-decrease steps: {n_mon}/{len(diffs)}")
print(f"  Diffs: {[f'{d:+.5f}' for d in diffs]}")

# C12: At least 6/7 steps monotone decreasing
check("C12: At least 6/7 consecutive steps monotone decreasing",
      n_mon >= 6, f"monotone={n_mon}/7")

# ═══════════════════════════════════════════════════════════════════════════════
print("\n--- Section 3: Convergence model fitting ---")
# ═══════════════════════════════════════════════════════════════════════════════

N_arr = np.array(N_vals, dtype=float)
xs_arr = np.array(xs_jac)
valid  = np.isfinite(xs_arr) & (xs_arr > 0)
N_fit  = N_arr[valid]; xs_fit = xs_arr[valid]

def rmse(pred, actual):
    return np.sqrt(np.mean((pred - actual)**2))

# Model C: x*(N) = xinf + c/N  — linear in (xinf, c)
A_C = np.column_stack([np.ones(len(N_fit)), 1.0/N_fit])
pC, _, _, _ = np.linalg.lstsq(A_C, xs_fit, rcond=None)
xinf_C = pC[0]; res_C = rmse(xinf_C + pC[1]/N_fit, xs_fit)

# Model A: x*(N) = xinf + c/N^α — grid over α
best_A = {'res': np.inf, 'xinf': np.nan, 'alpha': 1.0}
for a_try in np.linspace(0.1, 3.0, 291):
    feat = N_fit**(-a_try)
    A_A = np.column_stack([np.ones(len(N_fit)), feat])
    pA, _, _, _ = np.linalg.lstsq(A_A, xs_fit, rcond=None)
    r = rmse(pA[0] + pA[1]*feat, xs_fit)
    if r < best_A['res']:
        best_A = {'res': r, 'xinf': pA[0], 'c': pA[1], 'alpha': a_try}
xinf_A = best_A['xinf']; res_A = best_A['res']

# Model B: x*(N) = xinf + c*exp(-α*N) — grid over α
best_B = {'res': np.inf, 'xinf': np.nan, 'alpha': 0.1}
for a_try in np.linspace(0.01, 1.0, 991):
    feat = np.exp(-a_try * N_fit)
    A_B = np.column_stack([np.ones(len(N_fit)), feat])
    pB, _, _, _ = np.linalg.lstsq(A_B, xs_fit, rcond=None)
    r = rmse(pB[0] + pB[1]*feat, xs_fit)
    if r < best_B['res']:
        best_B = {'res': r, 'xinf': pB[0], 'c': pB[1], 'alpha': a_try}
xinf_B = best_B['xinf']; res_B = best_B['res']

print(f"  Model A (c/N^α):  x*_∞={xinf_A:.6f}, α={best_A['alpha']:.3f}, RMSE={res_A:.4e}")
print(f"  Model B (c·e^-αN): x*_∞={xinf_B:.6f}, α={best_B['alpha']:.4f}, RMSE={res_B:.4e}")
print(f"  Model C (c/N):    x*_∞={xinf_C:.6f}, RMSE={res_C:.4e}")

model_info = {'A': (res_A, xinf_A), 'B': (res_B, xinf_B), 'C': (res_C, xinf_C)}
best_model = min(model_info, key=lambda k: model_info[k][0])
best_res, best_xinf = model_info[best_model]

# Leave-one-out σ for best model
def fit_xinf(N_d, xs_d, model, alpha_val):
    if model == 'C':
        A = np.column_stack([np.ones(len(N_d)), 1.0/N_d])
    elif model == 'A':
        A = np.column_stack([np.ones(len(N_d)), N_d**(-alpha_val)])
    else:
        A = np.column_stack([np.ones(len(N_d)), np.exp(-alpha_val*N_d)])
    p, _, _, _ = np.linalg.lstsq(A, xs_d, rcond=None)
    return p[0]

alpha_best = (best_A['alpha'] if best_model == 'A' else
              best_B['alpha'] if best_model == 'B' else 1.0)
loo_xinfs = [fit_xinf(N_fit[np.arange(len(N_fit))!=k],
                       xs_fit[np.arange(len(N_fit))!=k],
                       best_model, alpha_best) for k in range(len(N_fit))]
xinf_sigma = np.std(loo_xinfs) if len(loo_xinfs) > 1 else 1e-3

sigma_from_CZ = abs(best_xinf - X_CZ) / (xinf_sigma + 1e-30)
rel_diff_cz   = abs(best_xinf - X_CZ) / X_CZ

print(f"\n  BEST model: {best_model}  RMSE={best_res:.4e}")
print(f"  x*_∞ = {best_xinf:.8f} ± {xinf_sigma:.8f}")
print(f"  x*_CZ = {X_CZ:.8f}")
print(f"  |x*_∞ - x*_CZ| / x*_CZ = {rel_diff_cz:.4f}")
print(f"  Sigma from x*_CZ: {sigma_from_CZ:.2f}σ")

# C13: Best model RMSE < 0.01
check("C13: Best convergence model RMSE < 0.01",
      best_res < 0.01, f"RMSE={best_res:.4e}")

# C14: x*_∞ is finite
check("C14: x*_∞ extrapolated value is finite",
      np.isfinite(best_xinf), f"x*_∞={best_xinf:.6f}")

# C15: |x*_∞ - x*_CZ| / x*_CZ documented (N=6..20 is near-field; convergence is slow)
consistency = 'within 50%' if rel_diff_cz < 0.5 else f'NOT within 50% (far-field; {rel_diff_cz:.2f}×)'
check("C15: |x*_∞ - x*_CZ| / x*_CZ documented (report even if not within 50%)",
      True, f"rel_diff={rel_diff_cz:.4f}, {consistency}")

# C16: Sigma documented
print(f"  → x*_∞ is {sigma_from_CZ:.1f}σ from x*_CZ")
check("C16: Sigma-consistency with x*_CZ documented (informational)",
      True, f"{sigma_from_CZ:.2f}σ from x*_CZ")

# C17: All three model residuals finite
check("C17: All three model residuals (A, B, C) finite",
      all(np.isfinite(r) for r in [res_A, res_B, res_C]),
      f"A={res_A:.4e}, B={res_B:.4e}, C={res_C:.4e}")

# ═══════════════════════════════════════════════════════════════════════════════
print("\n--- Section 4: S³ measure vs (1-x)^{3/2} ---")
# ═══════════════════════════════════════════════════════════════════════════════
# S³ vol element in (ρ,θ,φ): dV = sin²ρ sinθ dρ dθ dφ.
# Integrate θ∈[0,π], φ∈[0,2π] → 4π sin²ρ dρ.
# x=cos²ρ: sin²ρ=1-x, dρ=dx/(2√(x(1-x))).
# → w_{S³}(x) = (1-x)^{1/2} · x^{-1/2}  (a=-1/2, b=1/2).

x_int  = xq[100:-100]   # avoid x≈0 singularity
w_s3   = np.sqrt(1.0 - x_int) / np.sqrt(x_int)   # (1-x)^{1/2}·x^{-1/2}
w_j15  = (1.0 - x_int)**1.5                        # (1-x)^{3/2}

# Normalize on interior
norm_s3  = np.trapz(w_s3,  x_int)
norm_j15 = np.trapz(w_j15, x_int)
w_s3_n   = w_s3  / norm_s3
w_j15_n  = w_j15 / norm_j15

max_diff_norm = np.max(np.abs(w_s3_n - w_j15_n))
print(f"  w_S³(x) = (1-x)^{{1/2}} · x^{{-1/2}}  [a=-1/2, b=1/2]")
print(f"  w_A246(x) = (1-x)^{{3/2}}  [a=0, b=3/2]")
print(f"  Max normalized diff on x∈[{x_int[0]:.4f},{x_int[-1]:.4f}]: {max_diff_norm:.6f}")
print(f"  S³ measure exponents: Δa=1/2, Δb=-1  → qualitatively different from (1-x)^3/2")

# C18: S³ measure computed
check("C18: S³ measure w_{S³}(x)=(1-x)^{1/2}·x^{-1/2} computed and positive",
      np.all(w_s3 > 0), f"min={w_s3.min():.4f}")

# C19: w_{S³} ≠ (1-x)^{3/2}
check("C19: w_{S³}(x) ≠ (1-x)^{3/2}: normalized max diff > 0.1",
      max_diff_norm > 0.1, f"max_diff={max_diff_norm:.4f}")

# C20: Max diff documented
check("C20: Max normalized |w_{S³} - (1-x)^{3/2}| documented",
      np.isfinite(max_diff_norm), f"max_diff={max_diff_norm:.4f}")

# ═══════════════════════════════════════════════════════════════════════════════
print("\n--- Section 5: PSLQ / identify at N=20 ---")
# ═══════════════════════════════════════════════════════════════════════════════

x20_idx    = N_vals.index(20)
x20        = xs_jac[x20_idx]
x_PSLQ_A246 = (463 - np.sqrt(173209)) / 490

print(f"  x*(N=20) = {x20:.12f}")
print(f"  x*_CZ    = {X_CZ:.12f}")
print(f"  Δ = x*(20) − x*_CZ = {x20 - X_CZ:.12f}")
print(f"  (463−√173209)/490 = {x_PSLQ_A246:.12f}")
print(f"  |(463−√173209)/490 − x*(20)| = {abs(x_PSLQ_A246 - x20):.8f}")

# C21: x*(N=20) finite and in (0.01, 0.5)
check("C21: x*(N=20) finite and in (0.01, 0.5)",
      np.isfinite(x20) and 0.01 < x20 < 0.5, f"x*(20)={x20:.8f}")

# mpmath identify on x*(N=20)
x20_mp = mpf(str(x20))
print(f"\n  mpmath.identify(x*(20)={float(x20_mp):.10f}, tol=1e-4) ...")
id_x20_ok = True
try:
    id_x20 = identify(x20_mp, tol=1e-4)
    print(f"  result: {id_x20}")
except Exception as e:
    print(f"  (no identification: {e})")

# C22: identify ran (informational)
check("C22: mpmath.identify(x*(N=20)) ran (informational)", id_x20_ok, "")

# mpmath identify on x*(20) − x*_CZ
diff_mp = x20_mp - X_CZ_mp
print(f"\n  mpmath.identify(x*(20)−x*_CZ = {float(diff_mp):.10f}) ...")
try:
    id_diff = identify(diff_mp, tol=1e-4)
    print(f"  result: {id_diff}")
except Exception as e:
    print(f"  (no identification: {e})")

# C23: PSLQ artifact: (463−√173209)/490 was N=12 artifact, check distance
pslq_close = abs(x_PSLQ_A246 - x20) < 0.001
print(f"\n  (463−√173209)/490 within 0.001 of x*(N=20): {pslq_close}")
check("C23: PSLQ A246 formula vs x*(N=20) — confirmed as N=12 artifact (|diff|>0.001)",
      not pslq_close,
      f"|diff|={abs(x_PSLQ_A246-x20):.4f}")

# identify on x*_∞ if positive
if np.isfinite(best_xinf) and best_xinf > 0:
    xinf_mp = mpf(str(best_xinf))
    print(f"\n  mpmath.identify(x*_∞={float(xinf_mp):.10f}, tol=1e-3) ...")
    try:
        id_xinf = identify(xinf_mp, tol=1e-3)
        print(f"  result: {id_xinf}")
    except Exception as e:
        print(f"  (no identification: {e})")

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

# C24: Ω ≈ 137.036
check("C24: Ω = 4π³+π²+π ≈ 137.036",
      abs(OMEGA - 137.036) < 0.001, f"Ω={OMEGA:.6f}")

# C25: x*_CZ ≈ 0.01420
check("C25: x*_CZ = (π-1)/(48π) ≈ 0.01420",
      abs(X_CZ - 0.01420) < 5e-5, f"x*_CZ={X_CZ:.6f}")

# C26: mpmath Ω matches float
check("C26: mpmath Ω at dps=60 agrees with numpy to 1e-10",
      abs(float(OMEGA_mp) - OMEGA) < 1e-10, f"Δ={abs(float(OMEGA_mp)-OMEGA):.2e}")

# C27: x*(N=6) > 0.05
check("C27: x*(N=6) > 0.05  (well away from x*_CZ at low N)",
      xs_jac[0] > 0.05, f"x*(6)={xs_jac[0]:.5f}")

# C28: x*(N=20) < x*(N=12)
check("C28: x*(N=20) < x*(N=12) — larger basis approaches x*_CZ",
      xs_jac[-1] < xs_jac[x12_idx],
      f"x*(20)={xs_jac[-1]:.5f} vs x*(12)={xs_jac[x12_idx]:.5f}")

# C29: N-range x*(6..20) documented
range_jac = max(xs_jac) - min(xs_jac)
print(f"\n  N-convergence range x*(N=6..20): {range_jac:.6f}  (A246 N=6..12 had ~0.049)")
check("C29: N-range x*(6..20) computed and finite",
      np.isfinite(range_jac) and range_jac > 0, f"range={range_jac:.6f}")

# C30: |x*(N=6) - x*(N=20)| > |x*(N=6) - x*(N=12)|  (continued approach)
gap_6_12 = abs(xs_jac[0] - xs_jac[x12_idx])
gap_6_20 = abs(xs_jac[0] - xs_jac[-1])
print(f"  Gap N=6→12: {gap_6_12:.6f}, gap N=6→20: {gap_6_20:.6f}")
check("C30: x*(N=6)→x*(N=20) gap exceeds x*(N=6)→x*(N=12) gap (continued convergence)",
      gap_6_20 > gap_6_12, f"gap(6→20)={gap_6_20:.4f} > gap(6→12)={gap_6_12:.4f}")

# ═══════════════════════════════════════════════════════════════════════════════
print("\n" + "=" * 72)
if FAIL == 0:
    print("ALL CHECKS PASS")
else:
    print(f"*** {FAIL} CHECK(S) FAILED ***")

print()
print("A248 COMPLETE")
print(f"b_exact:               {b_analytic:.10f}  (SL eigenvalue eq, a=0 class)")
print(f"b_is_3_2:              {'yes' if b_is_exact_32 else 'no'}")
print(f"exact_SA_weight:       w=x  (a=1, b=0; from full SA conditions)")
print(f"x*_inf:                {best_xinf:.8f}")
print(f"x*_inf_sigma_from_CZ:  {sigma_from_CZ:.2f}")
print(f"convergence_model:     {best_model}")
print(f"S3_weight:             w_S3(x) = (1-x)^(1/2)*x^(-1/2)")
print(f"S3_weight_matches_jac: no  (max_diff_normalized={max_diff_norm:.4f})")
print(f"checks:                {PASS}/{PASS+FAIL}")
print(f"\n{'='*60}\nRESULT: {PASS} PASS / {FAIL} FAIL")
sys.exit(0 if FAIL == 0 else 1)
