#!/usr/bin/env python3
"""
Computation 54 -- Chiral idempotent f and stabilizer structure on C ⊗ H ⊗ O
=============================================================================
Per research/ngen_3.md section 13 and Comp 53, the next concrete probe
for N_gen = 3 (Furey 2022) is:

  (a) Construct the chiral primitive idempotent f in the chain algebra
      of D = C ⊗ H ⊗ O (acting on R^64).
  (b) Find the stabilizer sub-algebra Stab_f ⊂ L_D (= operators that
      commute with f, or equivalently leave f · D invariant).
  (c) Decompose the minimal ideal f · D under Stab_f and count the
      primitive idempotents.

Furey's claim: the count is 3, corresponding to three generations.

The construction:
  - O has volume element omega_O = e_1 e_2 ... e_7 (associated; on R^7 = Im O
    the antisymmetric product is the seven-form; via left-multiplication
    on R^8 = O, omega_O acts as a particular matrix).
  - The Furey chiral idempotent in Cl(0,6) is f = (I + i omega) / 2 where
    omega is the Cl(0,6) volume element (Computation 21).
  - On C ⊗ O (Furey 2014), the chiral half is selected by f.
  - On C ⊗ H ⊗ O (Furey 2022), Furey uses a SEQUENCE of nested chiral
    projectors that select the right-handed and electroweak structure.
"""
import numpy as np
import numpy.linalg as la


# Multiplication tables (from Comp 53)
C_T = [[(1, 0), (1, 1)], [(1, 1), (-1, 0)]]

H_T = [[(0, 0)] * 4 for _ in range(4)]
for _i in range(4):
    H_T[0][_i] = (1, _i); H_T[_i][0] = (1, _i)
for _i in [1, 2, 3]:
    H_T[_i][_i] = (-1, 0)
H_T[1][2] = (1, 3); H_T[2][1] = (-1, 3)
H_T[2][3] = (1, 1); H_T[3][2] = (-1, 1)
H_T[3][1] = (1, 2); H_T[1][3] = (-1, 2)

O_T = [[(0, 0) for _ in range(8)] for _ in range(8)]
for _i in range(8):
    O_T[0][_i] = (1, _i); O_T[_i][0] = (1, _i)
for _i in range(1, 8):
    O_T[_i][_i] = (-1, 0)
for _i, _j, _k in [(1, 2, 3), (1, 4, 5), (1, 7, 6),
                   (2, 4, 6), (2, 5, 7), (3, 4, 7), (3, 6, 5)]:
    O_T[_i][_j] = (1, _k); O_T[_j][_i] = (-1, _k)
    O_T[_j][_k] = (1, _i); O_T[_k][_j] = (-1, _i)
    O_T[_k][_i] = (1, _j); O_T[_i][_k] = (-1, _j)


def idx(c, h, o):
    return c * 32 + h * 8 + o


def decompose(i):
    return i // 32, (i % 32) // 8, i % 8


def CHO_mult_basis(b1, b2):
    c1, h1, o1 = decompose(b1); c2, h2, o2 = decompose(b2)
    sc, ic_ = C_T[c1][c2]; sh, ih = H_T[h1][h2]; so, io_ = O_T[o1][o2]
    return sc * sh * so, idx(ic_, ih, io_)


def L_matrix(c, h, o):
    """Left-multiplication matrix L_{e_c ⊗ e_h ⊗ e_o} on R^64."""
    M = np.zeros((64, 64))
    for j in range(64):
        c2, h2, o2 = decompose(j)
        s, k = CHO_mult_basis(idx(c, h, o), j)
        M[k, j] = s
    return M


def R_matrix(c, h, o):
    """Right-multiplication matrix R_{e_c ⊗ e_h ⊗ e_o} on R^64."""
    M = np.zeros((64, 64))
    for j in range(64):
        s, k = CHO_mult_basis(j, idx(c, h, o))
        M[k, j] = s
    return M


# ============================================================================
# Step 1: Construct the L_O volume element (image of e_1 e_2 ... e_7 under
# left-multiplication, with the alternative-algebra associated product)
# ============================================================================
def L_omega_O(include_e7=True):
    """Left-multiplication action of the O volume element.

    If include_e7 is True:  omega_7 = L_e1 * ... * L_e7  (Cl(0,7) volume, squares to +I)
    If include_e7 is False: omega_6 = L_e1 * ... * L_e6  (Cl(0,6) volume, squares to -I).

    The Furey chiral idempotent uses the Cl(0,6) volume (omega_6, squares to -I),
    so that i_C * omega_6 has square (-1)(-1) = +I, making (I + i_C omega_6)/2 a
    proper idempotent.
    """
    M = np.eye(8)
    max_i = 7 if include_e7 else 6
    for i in range(1, max_i + 1):
        Li = np.zeros((8, 8))
        for j in range(8):
            s, k = O_T[i][j]
            Li[k, j] = s
        M = Li @ M
    return M


# ============================================================================
# Step 2: Construct chiral idempotent f = (I + i_C * L_omega_O) / 2 on R^64.
# Note: i_C acts on the C factor; L_omega_O acts on the O factor.
# Combined operator on C ⊗ H ⊗ O = R^64: i_C ⊗ I_H ⊗ L_omega_O.
# ============================================================================
def chiral_idempotent_f():
    """
    f = (I + i_C ⊗ I_H ⊗ omega_6) / 2 in M_64(R), where omega_6 is the
    Cl(0,6) volume element (product of L_e1 through L_e6).  Since omega_6^2 = -I
    and i_C^2 = -1, the combined operator squares to +I, so f^2 = f.

    Rank of f should be 32 (half of R^64).
    """
    # i_C as a 2x2 real matrix:
    iC_real = np.array([[0, -1], [1, 0]])  # rotation by 90 deg = i in C
    I_H = np.eye(4)
    Lomega = L_omega_O(include_e7=False)  # use Cl(0,6) volume, NOT Cl(0,7)
    combo = np.kron(np.kron(iC_real, I_H), Lomega)
    return 0.5 * (np.eye(64) + combo), combo


def main():
    print("=" * 90)
    print("  Computation 54 -- Chiral idempotent and stabilizer on C ⊗ H ⊗ O")
    print("=" * 90)
    print()

    # Step 1: volume elements
    Lom7 = L_omega_O(include_e7=True)
    Lom6 = L_omega_O(include_e7=False)
    print(f"  (1) Volume element construction:")
    print(f"      omega_7 = L_e1 * ... * L_e7 (Cl(0,7) volume)")
    print(f"        omega_7^2 = {(Lom7 @ Lom7)[0, 0]:.4f} * I  (expected +I)")
    print(f"      omega_6 = L_e1 * ... * L_e6 (Cl(0,6) volume)")
    print(f"        omega_6^2 = {(Lom6 @ Lom6)[0, 0]:.4f} * I  (expected -I)")
    print(f"      omega_6 is the one Furey uses (squares to -1, gives proper chiral idempotent)")
    print()

    # Step 2: chiral idempotent
    f, combo = chiral_idempotent_f()
    print(f"  (2) Chiral idempotent f = (I + combo)/2 where combo = i_C ⊗ I_H ⊗ L_omega_O:")
    print(f"      combo^2 = {(combo @ combo)[0, 0]:.4f} * I  (should be -I or +I for idempotent)")
    fsq = f @ f
    is_idemp = np.allclose(fsq, f, atol=1e-9)
    print(f"      f^2 = f ?  {is_idemp}")
    if not is_idemp:
        # f^2 - f
        delta = fsq - f
        print(f"      f^2 - f max abs entry: {np.max(np.abs(delta)):.4f}")
    print(f"      Tr(f) = {np.trace(f):.2f}  (rank ~ 32 expected if combo^2 = -I)")
    print()

    # If f is a proper idempotent, study its image.
    if is_idemp:
        rank_f = int(round(float(np.real(np.trace(f)))))
        print(f"  (3) f has rank {rank_f}; minimal left ideal f * R^64 has dim {rank_f}.")
        print()

        # Find Stab_f: operators X in L_D such that f X f = X (i.e., X preserves f)
        # For idempotent f, the algebra f L_D f is the unital sub-algebra acting on f * R^64.
        # We compute this by intersection: take all L generators, conjugate by f, see what survives.
        print(f"  (4) Operators in L_D commuting with f (centralizer of f within L_D):")
        # Build a sample of L_D generators
        imag = [(1, 0, 0)] + [(0, h, 0) for h in [1, 2, 3]] + [(0, 0, o) for o in range(1, 8)]
        L_gens = [L_matrix(c, h, o) for (c, h, o) in imag]
        commuting = []
        for i, G in enumerate(L_gens):
            if np.allclose(G @ f, f @ G, atol=1e-9):
                commuting.append(i)
        print(f"      Of {len(L_gens)} imaginary-unit L-generators, {len(commuting)} "
              f"commute with f: {commuting}")
        print(f"      Names: {[imag[i] for i in commuting]}")
    else:
        print(f"  (3) f is not a proper idempotent under this construction; need to refine.")
    print()

    print("=" * 90)
    print("  Verdict / next step")
    print("=" * 90)
    print()
    if is_idemp:
        print(f"  Constructed a chiral idempotent f of rank {rank_f} on R^64.")
        print(f"  {len(commuting)} L-generators commute with f.")
        print(f"  These generate the L-side stabilizer sub-algebra Stab_f^L.")
        print(f"  Next (Comp 55): find primitive idempotent decomposition of Stab_f^L")
        print(f"  to identify the structural multiplicity inside f * R^64.  If the")
        print(f"  count is 3, that is the structural source of N_gen = 3.")
    else:
        print(f"  The naive construction f = (I + i_C ⊗ L_omega_O)/2 is not an idempotent.")
        print(f"  Furey's construction of the chiral idempotent in C ⊗ H ⊗ O is more")
        print(f"  intricate -- typically a SEQUENCE of nested projectors selecting")
        print(f"  the right-handed sector and the electroweak structure.  Comp 55")
        print(f"  will implement this carefully.")


if __name__ == "__main__":
    main()