#!/usr/bin/env python3
"""
Computation 43 -- Dixon-algebra foundation: chain algebra of O on C tensor O
==============================================================================
Per research/ngen_3.md section 11/12, Attack Line A proper requires
extending H_F = C^8 to a larger structure carrying the chain algebra
that Furey 2018 / 2022 identifies with three generations.

This script sets up the FOUNDATIONS:
  (1) Explicit octonion multiplication on R^8 via the Fano-plane.
  (2) Construction of left-action algebra L_O and right-action algebra R_O
      on O itself.  The chain algebra is the algebra of compositions
      L_a R_b for a, b in O.
  (3) Verification of L_O = R_O = Cl(0,7)_evenpart = Cl(0,6) on the
      complexification (Furey 2014 / Computation 16 result reproduced).
  (4) Probe of the right-action algebra R_O for natural 3-counts.
      Furey 2022 claims 3 generations emerge from the structure of
      mutually orthogonal primitive idempotents of R_O (or a related
      sub-algebra) acting on its own minimal ideal.

If a natural count of 3 emerges from R_O's idempotent structure, that is
the candidate structural source of N_gen = 3.  If not, this confirms
the structural conclusion that the count requires either input beyond
PST's postulates or further extension to the full Dixon algebra
R x C x H x O.
"""
import numpy as np
import numpy.linalg as la


# ============================================================================
# Step 1: Octonion multiplication via the Fano plane
# ============================================================================
# Basis: e_0 = 1 (real), e_1, ..., e_7 (imaginary units)
# Multiplication rules from the Fano plane (Cayley-Dickson convention):
#   e_i^2 = -1 for i = 1..7
#   e_i e_j = -e_j e_i for i != j (anti-commutative)
#   The 7 quaternionic triples (cyclic):
#     (1,2,3), (1,4,5), (1,7,6), (2,4,6), (2,5,7), (3,4,7), (3,6,5)
#   where (i,j,k) means e_i e_j = e_k, with cyclic permutations
#   negation under reverse: e_j e_i = -e_k

# Build the multiplication table
def build_octonion_table():
    """Return the structure constants C[i,j] = octonion product e_i * e_j
    encoded as (sign, index) where the product is sign * e_index."""
    n = 8
    table = [[(0, 0) for _ in range(n)] for _ in range(n)]
    # Real unit
    for i in range(n):
        table[0][i] = (1, i)
        table[i][0] = (1, i)
    # Diagonal: e_i^2 = -1 for i in 1..7
    for i in range(1, n):
        table[i][i] = (-1, 0)
    # Quaternionic triples (Fano-plane lines)
    triples = [
        (1, 2, 3),
        (1, 4, 5),
        (1, 7, 6),
        (2, 4, 6),
        (2, 5, 7),
        (3, 4, 7),
        (3, 6, 5),
    ]
    for i, j, k in triples:
        # e_i * e_j = +e_k, e_j * e_k = +e_i, e_k * e_i = +e_j
        # and the reverses are negative
        table[i][j] = (1, k); table[j][i] = (-1, k)
        table[j][k] = (1, i); table[k][j] = (-1, i)
        table[k][i] = (1, j); table[i][k] = (-1, j)
    return table


OCT_TABLE = build_octonion_table()


def octonion_mult(a, b):
    """Multiply two octonions given as length-8 real vectors."""
    out = np.zeros(8)
    for i in range(8):
        if abs(a[i]) < 1e-15:
            continue
        for j in range(8):
            if abs(b[j]) < 1e-15:
                continue
            sign, idx = OCT_TABLE[i][j]
            out[idx] += sign * a[i] * b[j]
    return out


def left_mult_matrix(a):
    """8x8 real matrix L_a: x -> a * x for x in O."""
    L = np.zeros((8, 8))
    for j in range(8):
        e_j = np.zeros(8); e_j[j] = 1.0
        L[:, j] = octonion_mult(a, e_j)
    return L


def right_mult_matrix(a):
    """8x8 real matrix R_a: x -> x * a for x in O."""
    R = np.zeros((8, 8))
    for j in range(8):
        e_j = np.zeros(8); e_j[j] = 1.0
        R[:, j] = octonion_mult(e_j, a)
    return R


def main():
    print("=" * 90)
    print("  Computation 43  --  Dixon-algebra foundation: chain algebra of O")
    print("=" * 90)
    print()

    # Sanity check: octonion multiplication
    print("  (1) Octonion multiplication, basic checks:")
    e = [np.eye(8)[i] for i in range(8)]
    # e_1 * e_1 = -1
    p = octonion_mult(e[1], e[1])
    print(f"    e_1 * e_1 = {p}  (expected [-1,0,...])")
    # e_1 * e_2 = e_3
    p = octonion_mult(e[1], e[2])
    print(f"    e_1 * e_2 = {p}  (expected [0,0,0,1,0,0,0,0])")
    # e_2 * e_1 = -e_3
    p = octonion_mult(e[2], e[1])
    print(f"    e_2 * e_1 = {p}  (expected [0,0,0,-1,0,0,0,0])")
    # Non-associativity: (e_1 e_2) e_4 vs e_1 (e_2 e_4)
    p1 = octonion_mult(octonion_mult(e[1], e[2]), e[4])
    p2 = octonion_mult(e[1], octonion_mult(e[2], e[4]))
    print(f"    (e_1 e_2) e_4 = {p1}")
    print(f"    e_1 (e_2 e_4) = {p2}")
    print(f"    Non-associative? {not np.allclose(p1, p2)}")
    print()

    # ============================================================================
    # Step 2: Left and right action algebras
    # ============================================================================
    print("  (2) Left- and right-action algebras L_O, R_O on R^8:")
    # Generate L_a for each basis e_i (i = 1..7; identity is trivial)
    L_basis = [left_mult_matrix(e[i]) for i in range(1, 8)]
    R_basis = [right_mult_matrix(e[i]) for i in range(1, 8)]
    # Check that L_O is associative (i.e., L_a L_b = L_{a*b}? NO -- L_O is associative
    # as MATRIX algebra but doesn't represent O faithfully because O is non-associative).
    # The right way: L_O is the algebra GENERATED by the L_{e_i} under matrix
    # multiplication.  It's known (Furey 2014) that this algebra is Cl(0,7).
    # Verify L_{e_i}^2 = -I for i = 1..7
    bad = 0
    for L in L_basis:
        if not np.allclose(L @ L, -np.eye(8)):
            bad += 1
    print(f"    L_{{e_i}}^2 = -I for i = 1..7: {7 - bad}/7 satisfied")
    # Verify anti-commutation L_{e_i} L_{e_j} + L_{e_j} L_{e_i} = -2 delta_{ij} I
    bad = 0
    for i in range(7):
        for j in range(i + 1, 7):
            ac = L_basis[i] @ L_basis[j] + L_basis[j] @ L_basis[i]
            if not np.allclose(ac, 0, atol=1e-10):
                bad += 1
    print(f"    Anti-commutation of L_{{e_i}}, L_{{e_j}} (i != j): {21 - bad}/21 satisfied")
    print(f"    => L_O contains 7 anti-commuting square-roots-of-(-I), so contains Cl(0,7)")
    print()

    # ============================================================================
    # Step 3: dimension of L_O as a matrix subalgebra
    # ============================================================================
    print("  (3) Dimension of L_O as a real matrix subalgebra of M_8(R):")
    # Generate all products of L_{e_i} up to length 7 and span them
    products = [np.eye(8)]
    for length in range(1, 5):  # Cl(0,7) acts on R^8 with at most dim 64; length 4 (7^4 = 2401 strings) saturates the span
        # Generate length-length products of basis L_{e_i}
        from itertools import product as it_product
        for indices in it_product(range(7), repeat=length):
            M = np.eye(8)
            for idx in indices:
                M = M @ L_basis[idx]
            products.append(M)
    # Flatten to vectors and find rank
    flat = np.stack([p.flatten() for p in products])
    rank_L = np.linalg.matrix_rank(flat, tol=1e-8)
    print(f"    rank of span(L_O products) in M_8(R) = {rank_L}")
    print(f"    (Cl(0,7) as real algebra has dim 128 = 2^7; rank should match if L_O ~= Cl(0,7))")
    print()

    # ============================================================================
    # Step 4: Probe right-action algebra R_O for natural 3-counts
    # ============================================================================
    print("  (4) Right-action algebra R_O: idempotent structure probe:")
    # Build complexified R_O acting on C^8 = R^8 + i R^8
    # Furey 2022's claim: R_O contains a copy of M_3(C) whose primitive idempotents
    # parameterize the three generations.
    # First check: R_O dim
    products_R = [np.eye(8)]
    for length in range(1, 5):  # Cl(0,7) acts on R^8 with at most dim 64; length 4 (7^4 = 2401 strings) saturates the span
        for indices in it_product(range(7), repeat=length):
            M = np.eye(8)
            for idx in indices:
                M = M @ R_basis[idx]
            products_R.append(M)
    flat_R = np.stack([p.flatten() for p in products_R])
    rank_R = np.linalg.matrix_rank(flat_R, tol=1e-8)
    print(f"    rank of span(R_O products) in M_8(R) = {rank_R}")
    print()

    # Find the center of R_O (operators commuting with all R_{e_i})
    # The center decomposes R_O into sectors; primitive idempotents of the center
    # count the "blocks" of R_O.
    print("    Center of R_O (operators in M_8(R) commuting with all R_{e_i}):")
    # Set up commutator-with-each-generator linear map
    rows = []
    for R in R_basis:
        # [R, M] in vec form: kron(R, I) - kron(I, R^T) applied to vec(M)
        L = np.kron(R, np.eye(8)) - np.kron(np.eye(8), R.T)
        rows.append(L)
    full = np.vstack(rows)
    null_dim = 64 - np.linalg.matrix_rank(full, tol=1e-8)
    print(f"    dim of center = {null_dim}")
    print()

    # Try the COMPLEXIFIED right-action algebra
    # On C^8 = C ⊗ R^8, R_O extends to a complex algebra; this is Furey's "complexified
    # chain algebra".  The primitive idempotent structure may reveal the 3-count.
    print("    Complexified R_O on C^8:")
    R_basis_C = [R.astype(complex) for R in R_basis]
    rows_C = []
    for R in R_basis_C:
        L = np.kron(R, np.eye(8, dtype=complex)) - np.kron(np.eye(8, dtype=complex), R.T)
        rows_C.append(L)
    full_C = np.vstack(rows_C)
    null_dim_C = 64 - np.linalg.matrix_rank(full_C, tol=1e-8)
    print(f"    dim of complex center = {null_dim_C}")
    print()

    print("=" * 90)
    print("  Verdict / next step")
    print("=" * 90)
    print()
    print("  The octonion left/right multiplication algebras are now constructed.")
    print(f"  Both L_O and R_O have matrix rank {rank_L} = {rank_R} = 64 = 2^{{log2(64)}} = dim M_8(R),")
    print("  so the chain algebra fills M_8(R) (consistent with Furey 2014 / Comp 16:")
    print("  O acts on itself via Cl(0,7) which is M_8(R) +/- M_8(R) at real level, M_8(C)")
    print("  at the complex level).")
    print()
    print("  For the 3-generation structure, Furey 2022 uses a SUBALGEBRA of R_O")
    print("  (not the full R_O) -- specifically, the operators that commute with a")
    print("  fixed primitive idempotent f (the Furey chiral projector).  The minimal")
    print("  ideal R_O · f under this subalgebra has the structure carrying the 3")
    print("  generations.")
    print()
    print("  Next concrete step (Comp 44 candidate): build the Furey primitive idempotent")
    print("  f in the C-O algebra setting, compute the stabilizer subalgebra Stab_f in R_O,")
    print("  and find the number of mutually orthogonal primitive idempotents in")
    print("  Stab_f acting on R_O · f.  If = 3, that is the Furey structural derivation")
    print("  of N_gen = 3 imported into PST.")


if __name__ == "__main__":
    main()