CCFU Script 2 — proof_g2_complete.py

"""
CCFU Proof 28 — Complete G₂(split) Identification
Exact rational arithmetic; no floating point computations.

The script computes the infinitesimal stabilizer g of the 3-form

    Ω = e0∧e1∧e4 + e0∧e2∧e5 + e0∧e3∧e6
        + e1∧e2∧e3 + e4∧e5∧e6

inside gl(7, R), and verifies the following statements exactly:

1. dim g = 14.
2. g is closed under commutators, hence a Lie subalgebra of gl(7, R).
3. [g, g] = g, so g is perfect.
4. The Killing form is nondegenerate, so g is semisimple by Cartan's criterion.
5. End_g(R^7) = R·Id. Since g is semisimple, complete reducibility implies
   that the 7-dimensional representation is irreducible. In fact this proves
   absolute irreducibility over R.
6. End_g(g_ad) = R·Id. For a semisimple real Lie algebra this implies that
   there is only one simple real ideal, hence g is simple.
7. The Killing form has exact signature (8, 6, 0). Using the classification of
   real simple Lie algebras / real forms of complex g₂, this identifies g as
   the split real form g₂(split), also denoted g_{2(2)}.
8. The bilinear form b_Ω computed directly from Ω has signature (3, 4, 0), and
   every X in g satisfies X^T b_Ω + b_Ω X = 0. Thus the computed algebra embeds
   in so(3,4), as expected for split g₂.

External mathematical inputs used in the final identification:
- Cartan criterion for semisimplicity.
- Complete reducibility of finite-dimensional representations of semisimple
  Lie algebras.
- Classification of real simple Lie algebras / real forms of complex g₂.

No appeal is made to Hitchin's classification theorem for stable 3-forms; the
bilinear form b_Ω is computed directly.

Run:
    python3 proof_g2_complete.py
"""

from fractions import Fraction
from itertools import combinations, permutations

N = 7
GL_DIM = N * N


def require(condition, message):
    """Runtime check that is not disabled by `python -O`."""
    if not condition:
        raise RuntimeError(message)


def sign_of_sequence(seq):
    """Sign of the permutation represented by `seq`, computed by inversions."""
    sign = 1
    items = list(seq)
    for i in range(len(items)):
        for j in range(i + 1, len(items)):
            if items[i] > items[j]:
                sign *= -1
    return sign


# === Build Ω ===
# Ω_W = τ∧ev + vol_W + vol_{W*} in the chosen ordered basis.
OMEGA_TERMS = [
    (0, 1, 4, +1),
    (0, 2, 5, +1),
    (0, 3, 6, +1),
    (1, 2, 3, +1),
    (4, 5, 6, +1),
]

Omega = {}
for i, j, k, coefficient in OMEGA_TERMS:
    for p in permutations((i, j, k)):
        Omega[p] = coefficient * sign_of_sequence(p)


def omega(i, j, k):
    return Omega.get((i, j, k), 0)


def rref(matrix, ncols):
    """Reduced row echelon form over Q, returning (rref_matrix, pivot_cols)."""
    M = [[Fraction(x) for x in row] for row in matrix]
    pivot_cols = []
    rank = 0

    for col in range(ncols):
        pivot = None
        for row in range(rank, len(M)):
            if M[row][col] != 0:
                pivot = row
                break
        if pivot is None:
            continue

        M[rank], M[pivot] = M[pivot], M[rank]
        pivot_value = M[rank][col]
        for c in range(col, ncols):
            M[rank][c] /= pivot_value

        for row in range(len(M)):
            if row != rank and M[row][col] != 0:
                factor = M[row][col]
                for c in range(col, ncols):
                    M[row][c] -= factor * M[rank][c]

        pivot_cols.append(col)
        rank += 1
        if rank == len(M):
            break

    return M, pivot_cols


def rank(matrix, ncols):
    return len(rref(matrix, ncols)[1])


def nullspace(matrix, ncols):
    """Basis of the kernel of `matrix`, using free variables as coordinates."""
    R, pivot_cols = rref(matrix, ncols)
    pivot_set = set(pivot_cols)
    free_cols = [c for c in range(ncols) if c not in pivot_set]

    basis = []
    for free_col in free_cols:
        v = [Fraction(0)] * ncols
        v[free_col] = Fraction(1)
        for row, pivot_col in enumerate(pivot_cols):
            v[pivot_col] = -R[row][free_col]
        basis.append(v)

    return basis, pivot_cols, free_cols


def v2m(v, size=N):
    return [[v[size * i + j] for j in range(size)] for i in range(size)]


def m2v(M):
    size = len(M)
    return [M[i][j] for i in range(size) for j in range(size)]


def matrix_mul(A, B):
    rows = len(A)
    inner = len(B)
    cols = len(B[0])
    return [
        [sum(A[i][k] * B[k][j] for k in range(inner)) for j in range(cols)]
        for i in range(rows)
    ]


def matrix_sub(A, B):
    return [
        [A[i][j] - B[i][j] for j in range(len(A[0]))]
        for i in range(len(A))
    ]


def commutator_vec(v1, v2):
    X = v2m(v1)
    Y = v2m(v2)
    return m2v(matrix_sub(matrix_mul(X, Y), matrix_mul(Y, X)))


def vector_linear_combination(coefficients, basis):
    if not basis:
        return []
    out = [Fraction(0)] * len(basis[0])
    for coeff, vector in zip(coefficients, basis):
        if coeff == 0:
            continue
        for i, value in enumerate(vector):
            out[i] += coeff * value
    return out


def exact_signature(M):
    """
    Exact inertia (pos, neg, zero) of a symmetric matrix over Q.

    This performs symmetric Gaussian elimination by congruence operations, so
    the inertia is preserved exactly by Sylvester's law of inertia.
    """
    n = len(M)
    S = [[Fraction(M[i][j]) for j in range(n)] for i in range(n)]
    diagonal_entries = []

    step = 0
    while step < n:
        pivot = None
        for i in range(step, n):
            if S[i][i] != 0:
                pivot = i
                break

        if pivot is None:
            # If all remaining diagonal entries are zero but an off-diagonal
            # entry is nonzero, add one row/column to another to create a
            # nonzero diagonal pivot. This is valid in characteristic != 2.
            for i in range(step, n):
                for j in range(i + 1, n):
                    if S[i][j] != 0:
                        for k in range(n):
                            S[i][k] += S[j][k]
                        for k in range(n):
                            S[k][i] += S[k][j]
                        pivot = i
                        break
                if pivot is not None:
                    break

        if pivot is None:
            break  # all remaining entries are zero

        if pivot != step:
            S[step], S[pivot] = S[pivot], S[step]
            for k in range(n):
                S[k][step], S[k][pivot] = S[k][pivot], S[k][step]

        d = S[step][step]
        require(d != 0, "internal error: zero pivot selected during signature computation")
        diagonal_entries.append(d)

        for i in range(step + 1, n):
            if S[i][step] == 0:
                continue
            factor = S[i][step] / d
            for j in range(step, n):
                S[i][j] -= factor * S[step][j]
            for j in range(step, n):
                S[j][i] -= factor * S[j][step]

        step += 1

    pos = sum(1 for d in diagonal_entries if d > 0)
    neg = sum(1 for d in diagonal_entries if d < 0)
    zero = n - len(diagonal_entries)
    return pos, neg, zero


def acts_trivially_on_omega(v):
    X = v2m(v)
    for i, j, k in BASIS_3:
        value = Fraction(0)
        for ell in range(N):
            value += X[ell][i] * omega(ell, j, k)
            value += X[ell][j] * omega(i, ell, k)
            value += X[ell][k] * omega(i, j, ell)
        if value != 0:
            return False
    return True


def commutant_dimension(matrices):
    """Dimension of {T : T A = A T for every A in matrices}."""
    size = len(matrices[0])
    equations = []
    for A in matrices:
        for row_index in range(size):
            for col_index in range(size):
                row = [Fraction(0)] * (size * size)
                for c in range(size):
                    # (T A)_{row_index,col_index}
                    row[size * row_index + c] += A[c][col_index]
                    # -(A T)_{row_index,col_index}
                    row[size * c + col_index] -= A[row_index][c]
                equations.append(row)
    return size * size - rank(equations, size * size)


def compute_b_omega():
    """
    Compute b_Ω exactly from

        b_Ω(u, v) vol = (1/6) (ι_u Ω) ∧ (ι_v Ω) ∧ Ω.

    The chosen volume form is e0∧...∧e6.
    """
    indices = list(range(N))
    b_omega = [[Fraction(0)] * N for _ in range(N)]

    for a in range(N):
        for b in range(N):
            total = Fraction(0)
            for c2a in combinations(indices, 2):
                remaining_after_a = [x for x in indices if x not in c2a]
                for c2b in combinations(remaining_after_a, 2):
                    c3 = tuple(x for x in remaining_after_a if x not in c2b)
                    if len(c3) != 3:
                        continue

                    va = omega(a, c2a[0], c2a[1])
                    vb = omega(b, c2b[0], c2b[1])
                    vo = omega(c3[0], c3[1], c3[2])
                    if va == 0 or vb == 0 or vo == 0:
                        continue

                    wedge_indices = list(c2a) + list(c2b) + list(c3)
                    if len(set(wedge_indices)) != N:
                        continue

                    total += sign_of_sequence(wedge_indices) * va * vb * vo

            b_omega[a][b] = total / Fraction(6)

    return b_omega


def preserves_bilinear_form(v, eta):
    X = v2m(v)
    for i in range(N):
        for j in range(N):
            value = sum(X[k][i] * eta[k][j] for k in range(N))
            value += sum(eta[i][k] * X[k][j] for k in range(N))
            if value != 0:
                return False
    return True


# === Step 1: infinitesimal stabilizer of Ω ===
BASIS_3 = [
    (i, j, k)
    for i in range(N)
    for j in range(i + 1, N)
    for k in range(j + 1, N)
]
GL_POSITIONS = [(a, b) for a in range(N) for b in range(N)]

stabilizer_equations = [[Fraction(0)] * GL_DIM for _ in range(len(BASIS_3))]
for row_index, (i, j, k) in enumerate(BASIS_3):
    for col_index, (a, b) in enumerate(GL_POSITIONS):
        value = 0
        if b == i:
            value += omega(a, j, k)
        if b == j:
            value += omega(i, a, k)
        if b == k:
            value += omega(i, j, a)
        stabilizer_equations[row_index][col_index] = Fraction(value)

basis_g, pivot_cols, free_cols = nullspace(stabilizer_equations, GL_DIM)
dim_g = len(basis_g)
print(f"Step 1: dim Stab(Ω) = {dim_g}")
require(dim_g == 14, f"Expected dim Stab(Ω) = 14, got {dim_g}")


def coordinates_in_g(v):
    """Coordinates in the stabilizer basis produced by nullspace()."""
    coefficients = [v[col] for col in free_cols]
    reconstructed = vector_linear_combination(coefficients, basis_g)
    require(reconstructed == v, "vector does not reconstruct from stabilizer coordinates")
    return coefficients


# === Step 2: closure under commutators ===
closed = all(
    acts_trivially_on_omega(commutator_vec(basis_g[i], basis_g[j]))
    for i in range(dim_g)
    for j in range(i + 1, dim_g)
)
print(f"Step 2: closed under commutators: {closed}")
require(closed, "The stabilizer is not closed under commutators")

# === Step 3: perfectness ===
derived_vectors = [
    commutator_vec(basis_g[i], basis_g[j])
    for i in range(dim_g)
    for j in range(i + 1, dim_g)
]
derived_dim = rank(derived_vectors, GL_DIM)
print(f"Step 3: dim [g,g] = {derived_dim} {'(PERFECT)' if derived_dim == dim_g else ''}")
require(derived_dim == dim_g, f"Expected [g,g] to have dimension {dim_g}, got {derived_dim}")

# === Structure constants in the computed basis of g ===
structure = [
    [[Fraction(0)] * dim_g for _ in range(dim_g)]
    for _ in range(dim_g)
]
for i in range(dim_g):
    for j in range(dim_g):
        if i == j:
            continue
        coeffs = coordinates_in_g(commutator_vec(basis_g[i], basis_g[j]))
        for k, coeff in enumerate(coeffs):
            structure[i][j][k] = coeff

# === Step 4: Killing form and semisimplicity ===
killing = [
    [
        sum(
            structure[i][a][b] * structure[j][b][a]
            for a in range(dim_g)
            for b in range(dim_g)
        )
        for j in range(dim_g)
    ]
    for i in range(dim_g)
]
killing_rank = rank(killing, dim_g)
print(f"Step 4: rank(Killing) = {killing_rank} {'(SEMISIMPLE)' if killing_rank == dim_g else ''}")
require(killing_rank == dim_g, "Killing form is degenerate; Cartan criterion does not give semisimplicity")

# === Step 5: natural commutant and irreducibility ===
natural_matrices = [v2m(v) for v in basis_g]
natural_commutant_dim = commutant_dimension(natural_matrices)
print(
    f"Step 5: dim End_g(R^7) = {natural_commutant_dim} "
    f"{'(ABSOLUTELY IRREDUCIBLE)' if natural_commutant_dim == 1 else ''}"
)
require(
    natural_commutant_dim == 1,
    f"Expected natural commutant dimension 1, got {natural_commutant_dim}",
)

# === Step 6: ad-commutant and simplicity ===
# ad(e_i) has entries row=target k, column=source j: [e_i,e_j] = sum_k c^k_ij e_k.
ad_matrices = []
for i in range(dim_g):
    ad_i = [[Fraction(0)] * dim_g for _ in range(dim_g)]
    for source in range(dim_g):
        for target in range(dim_g):
            ad_i[target][source] = structure[i][source][target]
    ad_matrices.append(ad_i)

ad_commutant_dim = commutant_dimension(ad_matrices)
print(f"Step 6: dim End_g(g_ad) = {ad_commutant_dim} {'(SIMPLE)' if ad_commutant_dim == 1 else ''}")
require(ad_commutant_dim == 1, f"Expected ad-commutant dimension 1, got {ad_commutant_dim}")

# === Step 7: exact Killing signature and split real form ===
k_pos, k_neg, k_zero = exact_signature(killing)
is_split_g2_signature = (k_pos, k_neg, k_zero) == (8, 6, 0)
print(
    f"Step 7: sig(Killing) = ({k_pos},{k_neg}), zero={k_zero} "
    f"{'(SPLIT G₂ SIGNATURE)' if is_split_g2_signature else ''}"
)
require(
    is_split_g2_signature,
    f"Expected exact split-g₂ Killing signature (8,6,0), got ({k_pos},{k_neg},{k_zero})",
)

# === Step 8: directly computed b_Ω and the so(3,4) embedding ===
eta = compute_b_omega()
metric_symmetric = all(eta[i][j] == eta[j][i] for i in range(N) for j in range(N))
require(metric_symmetric, "b_Ω is not symmetric")

eta_pos, eta_neg, eta_zero = exact_signature(eta)
print(f"Step 8a: sig(b_Ω) = ({eta_pos},{eta_neg}), zero={eta_zero}")
require(
    eta_zero == 0 and sorted((eta_pos, eta_neg)) == [3, 4],
    f"Expected b_Ω to have split signature (3,4,0) up to sign, got ({eta_pos},{eta_neg},{eta_zero})",
)

in_so_eta = all(preserves_bilinear_form(v, eta) for v in basis_g)
print(f"Step 8b: g ⊂ so(b_Ω): {in_so_eta} {'(so(3,4)-COMPATIBLE)' if in_so_eta else ''}")
require(in_so_eta, "g is not contained in so(b_Ω)")

print()
print("=" * 68)
print("CONCLUSION")
print("=" * 68)
print()
print("  Computed exact invariants:")
print("  1. dim g = 14")
print("  2. g is a Lie subalgebra of gl(7, R)")
print("  3. g is perfect: [g,g] = g")
print("  4. g is semisimple: Killing form is nondegenerate")
print("  5. R^7 is absolutely irreducible as a g-module")
print("  6. g is simple: the ad-commutant has dimension 1")
print("  7. sig(Killing) = (8,6), zero=0")
print("  8. b_Ω has signature (3,4), and g ⊂ so(b_Ω)")
print()
print("  Therefore, using the standard classification of real simple Lie")
print("  algebras / real forms of complex g₂, the computed Lie algebra is")
print("  isomorphic to the split real form g₂(split) = g_{2(2)}.")
print()
print("  VERIFIED with exact rational arithmetic.")
print("=" * 68)