#!/usr/bin/env python3
"""verify_P314.py -- Verifier for Addendum 314 (P18-T2 sub-result (i),
boundary Laplacian on S^3 + bulk Dirac^2 Lichnerowicz constant).

Recomputes from scratch and asserts, by machine, the two items A312 left
ARGUED rather than machine-checked:

  (A) The boundary Laplacian Delta_S^3 is the unique 2nd-order
      SO(4)-invariant operator up to scale. The Laplace-Beltrami eigenvalue
      on the degree-L harmonic block is -L(L+2), multiplicity (L+1)^2. The
      quadratic Casimir of so(4) = su(2)_L (+) su(2)_R, acting on the
      symmetric (j,j) irreps (L = 2j), equals 4 j(j+1) = L(L+2): it IS the
      Laplace-Beltrami operator up to sign/scale. The space of degree-<=2
      SO(4)-invariant operators is 2-dimensional, spanned by {I, Casimir}, so
      the only nonconstant 2nd-order invariant is the Laplacian up to affine
      a*Lap + b*I.

  (B) The bulk Dirac^2 Lichnerowicz CONSTANT R/4 = 3/2 on the round S^3.
      Lichnerowicz: D^2 = nabla*nabla + R/4 (scalar curvature R; R = 6 on the
      unit round S^3 -> R/4 = 3/2). The Dirac spectrum is +-(k+3/2), so D^2
      eigenvalues are (k+3/2)^2; the connection Laplacian has eigenvalues
      (k+3/2)^2 - 3/2; the difference is 3/2 = R/4, k-independent.

Honest scope: this machine-verifies the boundary Laplacian uniqueness UP TO
SCALE and the Dirac^2 Lichnerowicz CONSTANT. It does NOT machine-verify the
full bulk B^4 Dirac^2 operator (only the round-S^3 spectral identity) nor fix
the SCALE of the Laplacians. Sub-result (i) is ADVANCED, not closed; A294's
re-typing of the whole P18-T2 stands; (ii)=A310 chain and (iii)=A302 are
unchanged.

  S1  so(4) Casimir = S^3 Laplacian spectrum on the (j,j) irreps  - checks 1-4
  S2  Invariant 2nd-order operator space is 2-dim {I, Casimir}     - check 5
  S3  Dirac^2 Lichnerowicz constant R/4 = 3/2 on round S^3         - checks 6-8

Copyright Léon Fernando Vlegels -- CC BY 4.0
"""
import sys
import math

import numpy as np

PI = math.pi
PASS = FAIL = 0


def check(n, desc, cond):
    global PASS, FAIL
    ok = bool(cond)
    PASS += ok
    FAIL += (not ok)
    print(f"  [{'PASS' if ok else 'FAIL'}] {n:>2}. {desc}")


# --- su(2) spin-j generators from standard ladder operators ----------------
def su2(j):
    dim = int(round(2 * j + 1))
    ms = [j - i for i in range(dim)]
    Jp = np.zeros((dim, dim), dtype=complex)
    Jz = np.zeros((dim, dim), dtype=complex)
    for i, m in enumerate(ms):
        Jz[i, i] = m
        if i - 1 >= 0:
            Jp[i - 1, i] = math.sqrt(j * (j + 1) - m * (m + 1))
    Jm = Jp.conj().T
    return (Jp + Jm) / 2.0, (Jp - Jm) / (2.0j), Jz


def so4_gens(j):
    """[J_L^x,J_L^y,J_L^z, J_R^x,J_R^y,J_R^z] on the (j,j) tensor space, plus
    the identity."""
    Lx, Ly, Lz = su2(j)
    d = Lx.shape[0]
    I = np.eye(d, dtype=complex)
    g = [np.kron(A, I) for A in (Lx, Ly, Lz)]
    g += [np.kron(I, A) for A in (Lx, Ly, Lz)]
    return g, np.eye(d * d, dtype=complex)


def so4_casimir(g):
    """so(4) Casimir C2 = sum_i (A_i^2 + B_i^2), A_i=J_L+J_R, B_i=J_L-J_R
    = 2(J_L^2 + J_R^2). On (j,j): 4 j(j+1) = L(L+2), L=2j."""
    JL, JR = g[0:3], g[3:6]
    C = np.zeros_like(g[0])
    for i in range(3):
        A = JL[i] + JR[i]
        B = JL[i] - JR[i]
        C = C + A @ A + B @ B
    return C


# ===========================================================================
print("S1  so(4) Casimir on (j,j) equals the S^3 Laplacian spectrum L(L+2)")

# (1) su(2) Casimir J^2 = j(j+1) for each j
su2_ok = True
for j in (0.0, 0.5, 1.0, 1.5):
    Jx, Jy, Jz = su2(j)
    J2 = Jx @ Jx + Jy @ Jy + Jz @ Jz
    ev = np.linalg.eigvalsh((J2 + J2.conj().T) / 2.0).real
    if abs(float(np.mean(ev)) - j * (j + 1)) > 1e-9 or (ev.max() - ev.min()) > 1e-9:
        su2_ok = False
check(1, "su(2) Casimir J^2 = j(j+1) on each spin-j rep (j=0,1/2,1,3/2)",
      su2_ok)

# (2) so(4) Casimir on (j,j) is a scalar 4 j(j+1) (degenerate over the block)
scalar_ok = True
c2_vals = {}
for j in (0.0, 0.5, 1.0, 1.5):
    g, _ = so4_gens(j)
    C2 = so4_casimir(g)
    ev = np.linalg.eigvalsh((C2 + C2.conj().T) / 2.0).real
    c2_vals[j] = float(np.mean(ev))
    if (ev.max() - ev.min()) > 1e-8 or abs(c2_vals[j] - 4 * j * (j + 1)) > 1e-9:
        scalar_ok = False
check(2, "so(4) Casimir on (j,j) is the scalar 4 j(j+1) (block-degenerate)",
      scalar_ok)

# (3) 4 j(j+1) = L(L+2) with L = 2j -- equals the -Delta_S^3 spectrum
lap_ok = True
for j in (0.0, 0.5, 1.0, 1.5):
    L = int(round(2 * j))
    if abs(c2_vals[j] - L * (L + 2)) > 1e-9:
        lap_ok = False
check(3, "Casimir 4j(j+1) = L(L+2), L=2j: the -Delta_S^3 Laplacian spectrum "
      "(L=0->0, 1->3, 2->8, 3->15)", lap_ok)

# (4) block multiplicity (2j+1)^2 = (L+1)^2 -- the harmonic multiplicity
mult_ok = True
for j in (0.0, 0.5, 1.0, 1.5):
    g, I = so4_gens(j)
    L = int(round(2 * j))
    if I.shape[0] != (L + 1) ** 2:
        mult_ok = False
check(4, "degree-L harmonic multiplicity (2j+1)^2 = (L+1)^2 (1,4,9,16)",
      mult_ok)

# ===========================================================================
print("S2  SO(4)-invariant 2nd-order operator space is 2-dimensional {I, C2}")

# Build the j=1/2 (+) j=1 reducible space; on it I and C2 are independent
# (C2 = 3 on j=1/2, 8 on j=1). Enumerate degree-<=2 generator monomials,
# find the coefficient combinations that commute with all six generators
# (the invariant operators), and confirm the resulting operator span is
# exactly 2-dimensional = span{I, C2}.
blocks = [0.5, 1.0]
gen_list, I_list, C2_list, dims = [], [], [], []
for j in blocks:
    g, I = so4_gens(j)
    gen_list.append(g); I_list.append(I); C2_list.append(so4_casimir(g))
    dims.append(I.shape[0])
Dtot = sum(dims)


def embed(mats):
    M = np.zeros((Dtot, Dtot), dtype=complex)
    off = 0
    for A, d in zip(mats, dims):
        M[off:off + d, off:off + d] = A
        off += d
    return M


Iemb = embed(I_list)
C2emb = embed(C2_list)
Gemb = [embed([gen_list[b][a] for b in range(len(blocks))]) for a in range(6)]

monomials = [Iemb]
for a in range(6):
    monomials.append(Gemb[a])
for a in range(6):
    for b in range(6):
        monomials.append(Gemb[a] @ Gemb[b])
n_mono = len(monomials)

# coefficient combinations x with sum_i x_i [M_i, G_a] = 0 for all a
big = np.zeros((6 * Dtot * Dtot, n_mono), dtype=complex)
for i, m in enumerate(monomials):
    stacked = []
    for a in range(6):
        Ga = Gemb[a]
        stacked.append((m @ Ga - Ga @ m).reshape(-1))
    big[:, i] = np.concatenate(stacked)

u, s, vh = np.linalg.svd(big)
tol = 1e-7 * max(1.0, s[0]) if s.size else 1e-7
null_cols = np.where(s < tol)[0]
null_basis = vh.conj().T[:, null_cols] if s.size else np.zeros((n_mono, 0))
inv_ops = []
for c in range(null_basis.shape[1]):
    coeffs = null_basis[:, c]
    X = np.zeros((Dtot, Dtot), dtype=complex)
    for i in range(n_mono):
        X = X + coeffs[i] * monomials[i]
    inv_ops.append(X.reshape(-1))
inv_space_dim = (int(np.linalg.matrix_rank(np.array(inv_ops), tol=1e-7))
                 if inv_ops else 0)
ic2_rank = int(np.linalg.matrix_rank(
    np.array([Iemb.reshape(-1), C2emb.reshape(-1)]), tol=1e-8))

check(5, "invariant degree-<=2 operator space is 2-dim, spanned by {I, C2} "
      "(dim=%d, I&C2 independent rank=%d): Laplacian unique up to a*Lap+b*I"
      % (inv_space_dim, ic2_rank),
      inv_space_dim == 2 and ic2_rank == 2)

# ===========================================================================
print("S3  Bulk Dirac^2 Lichnerowicz constant R/4 = 3/2 on the round S^3")

R_S3 = 6.0                     # scalar curvature of the unit round S^3
lich = R_S3 / 4.0              # 3/2

# (6) R/4 = 3/2 on the unit round S^3
check(6, "scalar curvature R(S^3)=6, Lichnerowicz constant R/4 = 3/2",
      abs(lich - 1.5) < 1e-12)

# (7) Dirac^2 minus connection Laplacian = 3/2, k-independent (k=0..5)
# Dirac spectrum +-(k+3/2): D^2 = (k+3/2)^2; connection Laplacian on spinors
# (k+3/2)^2 - 3/2 (closed form). The difference recovers R/4.
const_ok = True
for k in range(0, 6):
    d2 = (k + 1.5) ** 2
    conn = (k + 1.5) ** 2 - 1.5
    if abs((d2 - conn) - lich) > 1e-12:
        const_ok = False
check(7, "D^2 - nabla*nabla = (k+3/2)^2 - [(k+3/2)^2 - 3/2] = 3/2 = R/4, "
      "k-independent (k=0..5)", const_ok)

# (8) the recovered constant matches R/4 exactly and is independent of the
# Dirac mode index (consistency: the spinor connection Laplacian spectrum
# (k+3/2)^2 - 3/2 is nonnegative for all k, as a connection Laplacian must be)
conn_nonneg = all(((k + 1.5) ** 2 - 1.5) >= -1e-12 for k in range(0, 6))
check(8, "recovered Lichnerowicz constant = R/4 = 3/2 and the connection "
      "Laplacian spectrum (k+3/2)^2 - 3/2 is nonnegative for all k",
      abs((((0 + 1.5) ** 2) - ((0 + 1.5) ** 2 - 1.5)) - lich) < 1e-12
      and conn_nonneg)

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