"""
verify_P044.py — P44: G₂ Action on J₃(𝕆) Off-Diagonals

Central numerical claims:
  1. rank(Bryant-3-form constraint matrix C) = 7  →  dim(G₂) = 14
  2. dim(Stab_{G₂}(e₇)) = 14 - 6 = 8  (stabilizer is SU(3))
  3. dim((S⁶)³ / G₂) = 3×6 - 14 = 4  mixing parameters
  4. θ_C = π/(2(h_{G₂}+1)) = π/14   (Cabibbo angle from Coxeter number h=6)
  5. Killing form eigenvalues for stabilizer algebra = negative & uniform (simple Lie alg)
"""

import numpy as np
import mpmath

mpmath.mp.dps = 50

PASS = FAIL = N = 0

def report(name, ok, claimed, actual):
    global PASS, FAIL, N
    ok = bool(ok); N += 1; PASS += ok; FAIL += (not ok)
    print(f"  [{'PASS' if ok else 'FAIL'}] {N:>2}. {name}")
    if not ok:
        print(f"          claimed={claimed}  actual={actual}")

pi = mpmath.pi

# Bryant 3-form (0-indexed: e₁→0,...,e₇→6)
# φ = e^{123}+e^{145}+e^{167}+e^{246}-e^{257}-e^{347}-e^{356}  (1-indexed)
# → (0-indexed):
fano_signed = [
    ((0,1,2),+1), ((0,3,4),+1), ((0,5,6),+1),
    ((1,3,5),+1), ((1,4,6),-1), ((2,3,6),-1), ((2,4,5),-1)
]
phi = np.zeros((7,7,7))
for (a,b,c),sgn in fano_signed:
    phi[a,b,c]=sgn; phi[b,c,a]=sgn; phi[c,a,b]=sgn
    phi[b,a,c]=-sgn; phi[a,c,b]=-sgn; phi[c,b,a]=-sgn

idx_pairs = [(i,j) for i in range(7) for j in range(i+1,7)]   # 21 entries
triples = [(i,j,k) for i in range(7) for j in range(i+1,7) for k in range(j+1,7)]  # 35

def vec_to_mat(v):
    A = np.zeros((7,7))
    for k,(i,j) in enumerate(idx_pairs):
        A[i,j]=v[k]; A[j,i]=-v[k]
    return A

# Build constraint matrix via Lie derivative:
# (L_A phi)_{ijk} = sum_l [A[l,i]*phi[l,j,k] + A[l,j]*phi[i,l,k] + A[l,k]*phi[i,j,l]]
C = np.zeros((35,21))
for row,(i,j,k) in enumerate(triples):
    for col in range(21):
        A = vec_to_mat(np.eye(21)[col])
        C[row,col] = sum(A[l,i]*phi[l,j,k]+A[l,j]*phi[i,l,k]+A[l,k]*phi[i,j,l] for l in range(7))

rank_C = np.linalg.matrix_rank(C, tol=1e-10)
dim_g2 = 21 - rank_C

# --- Check 1 ---
ok1 = (rank_C == 7) and (dim_g2 == 14)
report("P44-1: rank(C)=7 → dim(G₂)=14", ok1,
       "rank=7, dim=14", f"rank={rank_C}, dim={dim_g2}")

# G₂ generators
_,_,Vt = np.linalg.svd(C)
null_space = Vt[rank_C:].T
generators = [vec_to_mat(null_space[:,k]) for k in range(dim_g2)]

# --- Check 2: Stab(e₇) dim = 8 ---
e7 = np.zeros(7); e7[6]=1.0
M = np.column_stack([T@e7 for T in generators])
rank_M = np.linalg.matrix_rank(M, tol=1e-10)
dim_stab = dim_g2 - rank_M
ok2 = (rank_M == 6) and (dim_stab == 8)
report("P44-2: dim(Stab_{G₂}(e₇))=8", ok2,
       "rank_M=6, dim_stab=8", f"rank_M={rank_M}, dim_stab={dim_stab}")

# --- Check 3 ---
ok3 = (3*6 - 14 == 4)
report("P44-3: dim((S⁶)³/G₂)=3×6-14=4", ok3, 4, 3*6-14)

# --- Check 4 ---
h_G2 = 6
theta_C_formula = float(mpmath.pi / (2*(h_G2+1)))
theta_C_direct  = float(mpmath.pi / 14)
ok4 = abs(theta_C_formula - theta_C_direct) < 1e-15
report("P44-4: θ_C=π/(2(h_{G₂}+1))=π/14  (h_{G₂}=6)", ok4,
       f"π/14≈{theta_C_direct:.8f}", f"formula≈{theta_C_formula:.8f}")

# --- Check 5: Killing form for stabilizer algebra ---
# M has shape (7, 14): rows=e7 image components, cols=generator index
# null space of M (as linear map ℝ^14→ℝ^7) lives in ℝ^14
_,s2,Vt2 = np.linalg.svd(M)   # M is (7,14): U(7,7), s(7,), Vt(14,14)
# null space: rows of Vt2 with negligible singular values (index >= rank_M)
# But SVD gives s of length min(7,14)=7. Rows rank_M..14 of Vt2 are the null space.
stab_null = Vt2[rank_M:].T   # (14, 8) — components in generator basis
stab_gens = [sum(stab_null[k,j]*generators[k] for k in range(14)) for j in range(8)]

# Structure constants of stab algebra
n = 8
f_stab = np.zeros((n,n,n))
basis_flat = np.column_stack([g.flatten() for g in stab_gens])
for i in range(n):
    for j in range(i+1,n):
        comm = stab_gens[i]@stab_gens[j] - stab_gens[j]@stab_gens[i]
        coeffs,_,_,_ = np.linalg.lstsq(basis_flat, comm.flatten(), rcond=None)
        f_stab[i,j,:] = coeffs; f_stab[j,i,:] = -coeffs

K = np.array([[sum(f_stab[i,k,l]*f_stab[j,l,k] for k in range(n) for l in range(n))
               for j in range(n)] for i in range(n)])
evals = np.linalg.eigvalsh(K)
spread = max(evals)-min(evals)
mean_ev = np.mean(evals)
ok5 = (spread < 0.5) and (mean_ev < -1.0)  # all uniform & negative (simple compact)
report("P44-5: Killing form uniform & negative for Stab(e₇) [simple]", ok5,
       "all ≈ equal < 0",
       f"min={min(evals):.3f} max={max(evals):.3f} mean={mean_ev:.3f}")

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