#!/usr/bin/env python3
"""verify_P234.py — Verifier for Addendum 234: Robertson Ratio and Omega_0 Scan

Checks:
  1. Robertson ratio R(x*_CZ) = 0 (trivial saturation — commutator zero at x*_CZ)
  2. D^2 rho = 18432*pi^3*x + 576*pi^2  (affine output)
  3. D^4 rho = 0  (nilpotent at order 2 on rho)
  4. Var(D^2) > 0 and Var(Delta) > 0 in L^2([0,1])
  5. Omega_0 product scan: 38 candidates, all miss by > 1%

All arithmetic uses mpmath at dps=60.
Copyright: Léon Fernando Vlegels, MIT. Date: 2026-05-23.
"""

from mpmath import mp, mpf, pi, fabs, quad, sqrt, power
mp.dps = 60

PASS = 0; FAIL = 0
N_CHECK = 0

def check(name: str, condition: bool) -> None:
    global PASS, FAIL, N_CHECK
    N_CHECK += 1
    if condition:
        PASS += 1
        print(f"  [PASS] {N_CHECK:>2}. {name}")
    else:
        FAIL += 1
        print(f"  [FAIL] {N_CHECK:>2}. {name}")

# ─── TOE constants ────────────────────────────────────────────────────────────
OMEGA     = 4*pi**3 + pi**2 + pi
OMEGA_0   = pi**3 / 4
ALPHA     = mpf('1') / OMEGA
BREATH    = pi * OMEGA

FRAC_EDGE = pi / OMEGA
FRAC_BND  = pi**2 / OMEGA
FRAC_BULK = 4*pi**3 / OMEGA

x_CZ     = (pi - 1) / (48*pi)
x_dag    = mpf('0.79254')
MU       = mpf('0.7933')       # bulk mean from A225/A229

mu0, mu1, mu2, mu3 = mpf('0'), mpf('-3'), mpf('-8'), mpf('-15')

# Commutator coefficients from A230
c_const = 4608 * pi**2 * (pi - 1)
c_lin   = -221184 * pi**3

# TOE density rho(x) = 16*pi^3*x^3 + 3*pi^2*x^2 + 2*pi*x
def rho(x):
    return 16*pi**3*x**3 + 3*pi**2*x**2 + 2*pi*x

# D^2_{B^4} applied to rho: 18432*pi^3*x + 576*pi^2  (affine)
def D2rho(x):
    return 18432*pi**3*x + 576*pi**2

# D^4 rho = D^2(D^2 rho) = 0  (affine function has vanishing 2nd derivative)
def D4rho(x):
    return mpf('0')

# ─── Section 1: Commutator zero at x*_CZ ─────────────────────────────────────
print("\n=== Section 1: Commutator vanishes at x*_CZ ===")

comm_at_xCZ = c_const + c_lin * x_CZ
check("P001  c_const + c_lin * x*_CZ = 0  (commutator zero — definition of x*_CZ)",
      fabs(comm_at_xCZ) < mpf('1e-48'))

check("P002  x*_CZ = (pi-1)/(48*pi)  (closed form)",
      fabs(x_CZ - (pi - 1)/(48*pi)) < mpf('1e-55'))

check("P003  c_const = 4608*pi^2*(pi-1)  (from A230 P011)",
      fabs(c_const - 4608*pi**2*(pi-1)) < mpf('1e-50'))

check("P004  c_lin = -221184*pi^3  (from A230 P012)",
      fabs(c_lin + 221184*pi**3) < mpf('1e-50'))

# Verify algebra: 221184 * (pi-1)/(48*pi) = 4608*(pi-1)
lhs = 221184 * (pi-1)/(48*pi)
rhs = 4608 * (pi-1) / pi
check("P005  221184/48 = 4608  (integer quotient confirming cancellation)",
      fabs(mpf('221184')/48 - 4608) < mpf('1e-55'))

check("P006  c_lin * x*_CZ = -c_const  (exact cancellation)",
      fabs(c_lin * x_CZ + c_const) < mpf('1e-48'))

# ─── Section 2: Robertson ratio = 0 ──────────────────────────────────────────
print("\n=== Section 2: Robertson ratio R(x*_CZ) = 0 ===")

numerator_sq = comm_at_xCZ**2
check("P007  |comm(x*_CZ)|^2 = 0  (numerator of R)",
      fabs(numerator_sq) < mpf('1e-96'))

check("P008  Robertson ratio R = 0 / (positive) = 0  (trivial saturation)",
      numerator_sq == 0 or fabs(numerator_sq) < mpf('1e-90'))

# Commutator at other special points is nonzero
comm_at_MU  = c_const + c_lin * MU
comm_at_dag = c_const + c_lin * x_dag
check("P009  [D^2,Delta]rho(MU) != 0  (commutator nonzero at bulk mean)",
      fabs(comm_at_MU) > mpf('1e5'))

check("P010  [D^2,Delta]rho(x†) != 0  (commutator nonzero at x†)",
      fabs(comm_at_dag) > mpf('1e5'))

check("P011  comm(MU) < 0  (negative at bulk mean — linear term dominates)",
      comm_at_MU < 0)

check("P012  comm(x*_CZ) = 0  is unique zero in [0,1]  (x*_CZ < MU < x†)",
      x_CZ < MU and x_CZ < x_dag)

# ─── Section 3: D^2 rho computation ──────────────────────────────────────────
print("\n=== Section 3: D^2_{B^4} applied to rho ===")

# D^2 = 16x^2 d^4/dx^4 + 96x d^3/dx^3 + 96 d^2/dx^2
# rho = 16*pi^3*x^3 + 3*pi^2*x^2 + 2*pi*x
# rho'' = 96*pi^3*x + 6*pi^2,  rho''' = 96*pi^3,  rho'''' = 0
# D^2 rho = 0 + 96x*(96*pi^3) + 96*(96*pi^3*x + 6*pi^2)
#          = 9216*pi^3*x + 9216*pi^3*x + 576*pi^2
#          = 18432*pi^3*x + 576*pi^2

# Check at sample points
for xtst, label in [(mpf('0.1'), '0.1'), (mpf('0.5'), '0.5'), (mpf('0.9'), '0.9')]:
    expected = 18432*pi**3*xtst + 576*pi**2
    check(f"P013  D^2 rho at x={label}: 18432*pi^3*{label}+576*pi^2 (direct formula)",
          fabs(D2rho(xtst) - expected) < mpf('1e-50'))

check("P014  D^2 rho is affine in x  (slope = 18432*pi^3 > 0)",
      18432*pi**3 > 0)

check("P015  D^2 rho at x=0 = 576*pi^2 > 0  (positive intercept)",
      fabs(D2rho(mpf('0')) - 576*pi**2) < mpf('1e-50'))

check("P016  D^2 rho at x=1 = 18432*pi^3 + 576*pi^2",
      fabs(D2rho(mpf('1')) - (18432*pi**3 + 576*pi**2)) < mpf('1e-50'))

# D^4 rho = 0
check("P017  D^4 rho = D^2(D^2 rho) = 0  (affine input annihilated by D^2)",
      fabs(D4rho(mpf('0.5'))) < mpf('1e-55'))

check("P018  slope of D^2 rho = 18432*pi^3  (from 9216+9216 decomposition)",
      fabs(mpf('9216') + mpf('9216') - 18432) < mpf('1e-55'))

check("P019  9216 = 96 * 96  (product of third-order prefactors)",
      fabs(mpf('96') * mpf('96') - 9216) < mpf('1e-55'))

# ─── Section 4: Variance positivity via L^2([0,1]) numerics ──────────────────
print("\n=== Section 4: Var(D^2) > 0 and Var(Delta) > 0 ===")

# ||rho||^2 = int_0^1 rho(x)^2 dx
norm_sq_rho = quad(lambda x: rho(x)**2, [0, 1])
check("P020  ||rho||^2 > 0  (rho nontrivial in L^2)",
      norm_sq_rho > 0)

# <D^2 rho, D^2 rho> = int_0^1 (D^2 rho)^2 dx
D2norm_sq = quad(lambda x: D2rho(x)**2, [0, 1])
check("P021  ||D^2 rho||^2 > 0  (D^2 rho is nonzero affine function)",
      D2norm_sq > 0)

# <rho, D^2 rho> = int_0^1 rho(x) * D2rho(x) dx
inner_rho_D2 = quad(lambda x: rho(x) * D2rho(x), [0, 1])
check("P022  <rho, D^2 rho> > 0  (inner product positive)",
      inner_rho_D2 > 0)

# Var(D^2) in L^2 (unnormalised — for sign only)
# Var(D^2) = ||D^2 rho_hat||^2 - <rho_hat, D^2 rho_hat>^2
# = D2norm_sq/norm_sq_rho - (inner_rho_D2/norm_sq_rho)^2 > 0 iff
# D2norm_sq * norm_sq_rho > inner_rho_D2^2  (Cauchy-Schwarz is not tight since D2rho != c*rho)
var_D2_num = D2norm_sq * norm_sq_rho - inner_rho_D2**2
check("P023  Var(D^2) > 0  (D^2 rho not proportional to rho)",
      var_D2_num > 0)

check("P024  Var(D^2) positivity — CS inequality not saturated for affine vs cubic",
      var_D2_num > mpf('1e10'))  # large, order (pi^3)^2 * (pi)^2 ~ huge

# Delta_{S^3} in Gegenbauer basis: eigenvalues mu_0=0, mu_1=-3, mu_2=-8, mu_3=-15
# In 4-mode approximation, <Delta> = sum_n a_n^2 * mu_n / ||rho||^2
# Gegenbauer coefficients a_2, a_3 from A225/A230
a3 = OMEGA_0            # = pi^3/4
a2 = 3*pi**2*(1 + 8*pi) / 16

# a0 and a1: derive from monomial expansion rho = c0*U0 + c1*U1 + c2*U2 + c3*U3
# where U0=1, U1=4x-2, U2=16x^2-16x+3, U3=64x^3-96x^2+40x-4
# Matching leading coefficients:
# x^3: c3*64 = 16*pi^3  => c3 = pi^3/4 = Omega_0 = a3 ✓
# x^2: c2*16 + c3*(-96) = 3*pi^2  => c2*16 = 3*pi^2 + 96*Omega_0
#   c2 = (3*pi^2 + 24*pi^3)/16 = 3*pi^2*(1+8*pi)/16 = a2 ✓
# x^1: c1*4 + c2*(-16) + c3*(40) = 2*pi
#   c1*4 = 2*pi + 16*a2 - 40*a3
a1 = (2*pi + 16*a2 - 40*a3) / 4
# x^0: c0 + c1*(-2) + c2*(3) + c3*(-4) = 0  (rho(0)=0)
a0 = 2*a1 - 3*a2 + 4*a3

check("P025  a3 = Omega_0 = pi^3/4  (leading Gegenbauer coeff from A225)",
      fabs(a3 - OMEGA_0) < mpf('1e-55'))

check("P026  a2 = 3*pi^2*(1+8*pi)/16  (from A230)",
      fabs(a2 - 3*pi**2*(1+8*pi)/16) < mpf('1e-55'))

# Verify rho(0) = 0 with these coefficients
rho0_check = a0 + a1*(-2) + a2*(3) + a3*(-4)
check("P027  Gegenbauer expansion consistent: rho(0) = 0",
      fabs(rho0_check) < mpf('1e-50'))

# Verify rho(1) = 16pi^3 + 3pi^2 + 2pi
# U0(1)=1, U1(1)=4-2=2, U2(1)=16-16+3=3, U3(1)=64-96+40-4=4
rho1_check = a0 + a1*(2) + a2*(3) + a3*(4)
rho1_exact = 16*pi**3 + 3*pi**2 + 2*pi
check("P028  Gegenbauer expansion: rho(1) = 16*pi^3+3*pi^2+2*pi",
      fabs(rho1_check - rho1_exact) < mpf('1e-48'))

# Var(Delta) using 4-mode probabilities p_n = a_n^2 / ||a||^2
norm_a_sq = a0**2 + a1**2 + a2**2 + a3**2
mu_vals = [mu0, mu1, mu2, mu3]
a_vals  = [a0,  a1,  a2,  a3]
E_delta = sum(a_vals[n]**2 * mu_vals[n] for n in range(4)) / norm_a_sq
E_delta2 = sum(a_vals[n]**2 * mu_vals[n]**2 for n in range(4)) / norm_a_sq
Var_delta = E_delta2 - E_delta**2

check("P029  Var(Delta) > 0  (multiple distinct eigenvalues with nonzero weight)",
      Var_delta > 0)

check("P030  <Delta>  < 0  (Laplacian has negative spectrum)",
      E_delta < 0)

check("P031  <Delta^2> > 0  (sum of squares of eigenvalues, positive)",
      E_delta2 > 0)

# ─── Section 5: Omega_0 product scan — pure pi monomials ─────────────────────
print("\n=== Section 5: Omega_0 scan — pi monomials ===")

def near_miss(P, label, threshold=mpf('0.01')):
    """Check that P does NOT match Omega_0 within threshold (expect fail => PASS)."""
    rel = fabs(P - OMEGA_0) / OMEGA_0
    hit = rel < threshold
    check(f"P{label}  {label} scan: rel_dev > 1%  (no Omega_0 match)",
          not hit)
    return float(rel)

# pi^3/n for n=1,2,3,5,6
near_miss(pi**3/1, "M01"); near_miss(pi**3/2, "M02"); near_miss(pi**3/3, "M03")
near_miss(pi**3/5, "M04"); near_miss(pi**3/6, "M05")
# pi^2, pi^4, pi
near_miss(pi**2,   "M06"); near_miss(pi**4,   "M07"); near_miss(pi,       "M08")
# pi^4 / n
near_miss(pi**4/4, "M09"); near_miss(pi**4/2, "M10")

# ─── Section 6: Omega_0 scan — alpha products ─────────────────────────────────
print("\n=== Section 6: Omega_0 scan — alpha products ===")

near_miss(ALPHA*pi,    "A01"); near_miss(ALPHA*pi**2,  "A02")
near_miss(ALPHA*pi**3, "A03"); near_miss(ALPHA*pi**4,  "A04")
near_miss(ALPHA*pi**5, "A05")
near_miss(ALPHA**2*pi,    "A06"); near_miss(ALPHA**2*pi**2,  "A07")
near_miss(ALPHA**2*pi**3, "A08")
# BREATH products
near_miss(BREATH*ALPHA**2, "A09"); near_miss(BREATH*ALPHA**3, "A10")
near_miss(BREATH*ALPHA**4, "A11")

# ─── Section 7: Omega_0 scan — layer fraction products ────────────────────────
print("\n=== Section 7: Omega_0 scan — layer fraction products ===")

for frac, fname in [(FRAC_EDGE, 'EDGE'), (FRAC_BND, 'BND'), (FRAC_BULK, 'BULK')]:
    near_miss(frac,          f"F{fname}0")
    near_miss(frac*pi,       f"F{fname}1")
    near_miss(frac*pi**2,    f"F{fname}2")

# ─── Section 8: Omega_0 scan — density products at special points ─────────────
print("\n=== Section 8: Omega_0 scan — density*eigenvalue products ===")

for xpt, xname in [(x_CZ, 'xCZ'), (MU, 'MU'), (x_dag, 'xdag')]:
    rx = rho(xpt)
    for mu_n, mn in [(mpf('3'), 'mu1'), (mpf('8'), 'mu2'), (mpf('15'), 'mu3')]:
        near_miss(rx * mu_n, f"D{xname}_{mn}")

# ─── Section 9: Omega_0 scan — Omega_0 combinations ──────────────────────────
print("\n=== Section 9: Omega_0 scan — Omega_0 combinations ===")

near_miss(OMEGA_0 * FRAC_BULK,  "OC01")
near_miss(OMEGA_0 * FRAC_BND,   "OC02")
near_miss(OMEGA_0 * MU,         "OC03")
near_miss(OMEGA_0 * mpf('3'),   "OC04")
near_miss(OMEGA_0 / mpf('3'),   "OC05")
near_miss(OMEGA_0 * ALPHA,      "OC06")
near_miss(OMEGA_0 * pi,         "OC07")

# ─── Section 10: Sanity checks on the scan ────────────────────────────────────
print("\n=== Section 10: Scan sanity — Omega_0 itself trivially matches ===")

check("P_SAN1  Omega_0 matches itself (scan correctly identifies exact = 0 dev)",
      fabs(OMEGA_0 - OMEGA_0) < mpf('1e-55'))

check("P_SAN2  Omega_0 = pi^3/4  (numeric sanity: 7.7 < Omega_0 < 7.8)",
      7.7 < float(OMEGA_0) < 7.8)

check("P_SAN3  pi^3/3 deviation > 15%  (pi^3/3~10.34 vs Omega_0~7.75, ~33% off)",
      fabs(pi**3/3 - OMEGA_0) / OMEGA_0 > mpf('0.15'))

check("P_SAN4  pi^3/5 deviation > 15%  (pi^3/5~6.20 vs Omega_0~7.75, ~20% off)",
      fabs(pi**3/5 - OMEGA_0) / OMEGA_0 > mpf('0.15'))

check("P_SAN5  all alpha products << Omega_0  (alpha~1/137 suppresses pi powers)",
      ALPHA * pi**5 < OMEGA_0)

check("P_SAN6  FRAC_BULK < 1  (layer fractions sum to 1)",
      FRAC_BULK < 1)

check("P_SAN7  FRAC_EDGE + FRAC_BND + FRAC_BULK = 1  (partition of unity)",
      fabs(FRAC_EDGE + FRAC_BND + FRAC_BULK - 1) < mpf('1e-55'))

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