CCFU Script 1 — verify_higher_closure.py

id
2605225741400
title
CCFU Script 1 — verify_higher_closure.py
date
05/22/2026
text
Show source code
"""
verify_higher_closure.py
========================
Verifies the internal computational claims of Proof 21 (Higher Alternating Closure).
All computations exact over Q (Fraction). No numpy. No float.
No external dependencies — Poly class is implemented inline.

Conventions
-----------
DIM = 7 is the ambient dimension of V_W = R ⊕ W ⊕ W*.
Coordinates: index 0 = τ (R-summand), 1,2,3 = W, 4,5,6 = W*.
3-forms are stored as dicts {(i,j,k): coeff}, populated on all permutations
of every triple with appropriate signs, so a lookup by any tuple works.

Action convention for an X ∈ gl(V) on a 3-form ω:
  (X·ω)(i,j,k) = Σ_l X[l][i]·ω(l,j,k) + X[l][j]·ω(i,l,k) + X[l][k]·ω(i,j,l).
This is the standard Lie-derivative convention (free index = position in form,
summed index = first index of the matrix). action_on_lambda3,
infinitesimal_stabilizer_constraints, and the Part 6 inline check
all use this exact convention.

Parts
-----
  1. dim SL(W)-invariant 3-forms on V_W = R⊕W⊕W* equals 3.
  2. b_Ω/6 has the predicted block structure
     [symbolic proof in (a,b,c) over Q, plus numerical illustrations];
     det(b_Ω/6) is also checked directly symbolically, hence
     det(b_Ω) = 6⁷ · a⁹b⁶c⁶ / 64 = 4374·a⁹b⁶c⁶.
  3. sig(b_Ω) = (3,4) and dim Stab(Ω_{1,1,1}) = 14 — both exact.
     The identification with G₂_split uses the standard classification theorem
     for stable 3-forms in dimension 7.
  4. Ω_W and Ω_old as 3-forms on R⁷: exact integer witness H with
     H*Ω_W = 2·Ω_old; rescaling by 2^(-1/3) gives exact GL(7,R)
     equivalence. A determinant obstruction proves there is no exact
     GL(7,Q)-equivalence Ω_W ~ Ω_old.
  5. On W⊕W* (dim 6): dim Stab = 16 (exact). Not G₂.
  6. SL(W) preserves Ω_W: all 8 sl(3) generators verified exact.

Explicit runtime checks use require(), so they still run under python3 -O.
Pure Python + fractions.Fraction.

Run: python3 verify_higher_closure.py
"""

from fractions import Fraction
from itertools import permutations, combinations

DIM = 7          # ambient dimension of V_W
ZERO = Fraction(0)
ONE  = Fraction(1)


def require(condition, message):
    """Like assert, but unaffected by python3 -O (which strips assertions)."""
    if not condition:
        raise AssertionError(message)


# ============================================================
# Generic linear-algebra / combinatorial helpers
# ============================================================

def permutation_sign(seq):
    """Sign (parity of inversions) of a permutation given as a sequence."""
    items = list(seq)
    sign = 1
    length = len(items)
    for i in range(length):
        for j in range(i + 1, length):
            if items[i] > items[j]:
                sign = -sign
    return sign


def sorted_triples(dim):
    """All sorted triples (i,j,k) with 0 ≤ i < j < k < dim."""
    return [(i, j, k)
            for i in range(dim)
            for j in range(i + 1, dim)
            for k in range(j + 1, dim)]


def antisymmetric_form(signed_sorted):
    """
    Build a 3-form dict from values on sorted triples.

    signed_sorted: {(i,j,k): value} with i<j<k.
    `value` may be int / Fraction / Poly — anything that supports * by ±1.
    Raw ints are promoted to Fraction so downstream arithmetic stays exact.
    """
    form = {}
    for sorted_triple, value in signed_sorted.items():
        if isinstance(value, int):
            base = Fraction(value)
        else:
            base = value
        for perm in permutations(sorted_triple):
            form[perm] = base * permutation_sign(perm)
    return form


def zero_matrix(rows, cols=None):
    """rows × cols matrix of ZERO. Square if cols omitted."""
    if cols is None:
        cols = rows
    return [[ZERO] * cols for _ in range(rows)]


def row_reduce(matrix, num_cols):
    """Reduced row-echelon over Q (in place). Returns rank."""
    rank = 0
    num_rows = len(matrix)
    for col in range(num_cols):
        pivot_row = None
        for row in range(rank, num_rows):
            if matrix[row][col] != ZERO:
                pivot_row = row
                break
        if pivot_row is None:
            continue
        matrix[rank], matrix[pivot_row] = matrix[pivot_row], matrix[rank]
        for row in range(num_rows):
            if row != rank and matrix[row][col] != ZERO:
                factor = matrix[row][col] / matrix[rank][col]
                for c in range(num_cols):
                    matrix[row][c] -= factor * matrix[rank][c]
        rank += 1
    return rank


def kernel_dimension(stacked_rows, num_cols):
    """Dim of the kernel of a system Ax = 0 given by stacked rows."""
    work = [row[:] for row in stacked_rows]
    return num_cols - row_reduce(work, num_cols)


def exact_det(matrix, n):
    """Determinant by Leibniz over Q or any compatible exact commutative ring."""
    total = ZERO
    for perm in permutations(range(n)):
        sign = permutation_sign(perm)
        product = ONE
        for i in range(n):
            product *= matrix[i][perm[i]]
        total += sign * product
    return total


def prime_factorization_abs_integer(n):
    """Prime factorization of |n| as {prime: exponent}; n must be nonzero."""
    n = abs(int(n))
    if n == 0:
        raise ValueError("factorization of zero is undefined")
    factors = {}
    divisor = 2
    while divisor * divisor <= n:
        while n % divisor == 0:
            factors[divisor] = factors.get(divisor, 0) + 1
            n //= divisor
        divisor += 1 if divisor == 2 else 2   # after 2, test odd divisors only
    if n > 1:
        factors[n] = factors.get(n, 0) + 1
    return factors


def is_rational_nth_power(value, exponent):
    """
    True iff value ∈ Q is an exponent-th power of a rational number.

    For value = m/n in lowest terms, this is equivalent to every prime
    exponent in |m| and |n| being divisible by `exponent`, with the usual
    sign obstruction for even exponent.
    """
    if exponent <= 0:
        raise ValueError("exponent must be positive")

    value = Fraction(value)
    if value == 0:
        return True
    if value < 0 and exponent % 2 == 0:
        return False

    for factorization_target in (value.numerator, value.denominator):
        for power in prime_factorization_abs_integer(factorization_target).values():
            if power % exponent != 0:
                return False
    return True


def exact_signature(symmetric_matrix, n):
    """
    Signature (p, q) of a symmetric n×n matrix over Q, computed by
    exact congruence diagonalization (LDLᵀ with full pivoting).
    """
    S = [[symmetric_matrix[i][j] for j in range(n)] for i in range(n)]
    diagonals = []

    for step in range(n):
        pivot = -1
        for i in range(step, n):
            if S[i][i] != ZERO:
                pivot = i
                break

        if pivot == -1:
            found = False
            for i in range(step, n):
                for j in range(i + 1, n):
                    if S[i][j] != ZERO:
                        for k in range(n):
                            S[i][k] += S[j][k]
                        for k in range(n):
                            S[k][i] += S[k][j]
                        pivot = i
                        found = True
                        break
                if found:
                    break
            if not found:
                break

        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]

        diag = S[step][step]
        if diag == ZERO:
            continue
        diagonals.append(diag)

        for i in range(step + 1, n):
            if S[i][step] == ZERO:
                continue
            factor = S[i][step] / diag
            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]

    positive = sum(1 for d in diagonals if d > 0)
    negative = sum(1 for d in diagonals if d < 0)
    return positive, negative


# ============================================================
# Multivariate polynomial in (a, b, c) over Q
# ============================================================

class Poly:
    """
    Multivariate polynomial in three variables a, b, c over Q.

    Stored as a dict {(deg_a, deg_b, deg_c): Fraction}, with zero coefficients
    pruned. Supports +, -, *, **, /scalar, ==.
    Acts as a drop-in replacement for Fraction in the verifier's arithmetic,
    so build_omega_abc and hitchin_form_b_over_6 can run symbolically without
    any change.

    Two polynomials are equal iff their pruned dicts are equal — this is a
    rigorous polynomial-identity check, not a sample-based one.
    """
    __slots__ = ('terms',)

    def __init__(self, terms=None):
        if terms is None:
            self.terms = {}
        else:
            self.terms = {k: Fraction(v) for k, v in terms.items() if v != 0}

    @staticmethod
    def constant(value):
        v = Fraction(value)
        return Poly({(0, 0, 0): v}) if v != 0 else Poly()

    @staticmethod
    def variable(index):
        """0 = a, 1 = b, 2 = c."""
        exponents = [0, 0, 0]
        exponents[index] = 1
        return Poly({tuple(exponents): Fraction(1)})

    # ----- arithmetic -----

    def _coerce(self, other):
        if isinstance(other, Poly):
            return other
        return Poly.constant(other)

    def __add__(self, other):
        other = self._coerce(other)
        out = dict(self.terms)
        for k, v in other.terms.items():
            new = out.get(k, Fraction(0)) + v
            if new == 0:
                out.pop(k, None)
            else:
                out[k] = new
        result = Poly()
        result.terms = out
        return result

    def __radd__(self, other):
        return self.__add__(other)

    def __neg__(self):
        return Poly({k: -v for k, v in self.terms.items()})

    def __sub__(self, other):
        return self + (-self._coerce(other))

    def __rsub__(self, other):
        return self._coerce(other) - self

    def __mul__(self, other):
        other = self._coerce(other)
        out = {}
        for k1, v1 in self.terms.items():
            for k2, v2 in other.terms.items():
                k = (k1[0] + k2[0], k1[1] + k2[1], k1[2] + k2[2])
                new = out.get(k, Fraction(0)) + v1 * v2
                if new == 0:
                    out.pop(k, None)
                else:
                    out[k] = new
        result = Poly()
        result.terms = out
        return result

    def __rmul__(self, other):
        return self.__mul__(other)

    def __pow__(self, exponent):
        if not isinstance(exponent, int) or exponent < 0:
            raise ValueError("Poly power must be non-negative int")
        if exponent == 0:
            return Poly.constant(1)
        result = self
        for _ in range(exponent - 1):
            result = result * self
        return result

    def __truediv__(self, other):
        # Only scalar division. Division by a non-constant Poly would in
        # general leave Q[a,b,c], so we don't implement it.
        if isinstance(other, Poly):
            if len(other.terms) == 1 and (0, 0, 0) in other.terms:
                divisor = other.terms[(0, 0, 0)]
            else:
                raise ValueError("Poly division by non-constant unsupported")
        else:
            divisor = Fraction(other)
        if divisor == 0:
            raise ZeroDivisionError("division by zero polynomial")
        return Poly({k: v / divisor for k, v in self.terms.items()})

    # ----- comparison -----

    def __eq__(self, other):
        if not isinstance(other, Poly):
            other = Poly.constant(other)
        return self.terms == other.terms

    def __ne__(self, other):
        return not self.__eq__(other)

    def __bool__(self):
        return bool(self.terms)

    def __hash__(self):
        return hash(frozenset(self.terms.items()))

    def __repr__(self):
        if not self.terms:
            return "0"
        chunks = []
        for exponents in sorted(self.terms.keys(), reverse=True):
            coef = self.terms[exponents]
            mono_parts = []
            for var, deg in zip("abc", exponents):
                if deg == 1:
                    mono_parts.append(var)
                elif deg > 1:
                    mono_parts.append(f"{var}^{deg}")
            mono = "·".join(mono_parts)
            if not mono:
                chunks.append(f"{coef}")
            elif coef == 1:
                chunks.append(mono)
            elif coef == -1:
                chunks.append(f"-{mono}")
            else:
                chunks.append(f"{coef}·{mono}")
        return " + ".join(chunks).replace("+ -", "- ")


# ============================================================
# Lie-algebra structure: sl(3) and its lift to gl(V_W)
# ============================================================

def sl3_generators():
    """
    Eight standard generators of sl(3, R):
      six off-diagonal E_{ij} (i ≠ j),
      plus H₁ = E₁₁ - E₂₂, H₂ = E₂₂ - E₃₃.
    """
    gens = []
    for i in range(3):
        for j in range(3):
            if i != j:
                E = zero_matrix(3, 3)
                E[i][j] = ONE
                gens.append(E)

    H1 = zero_matrix(3, 3); H1[0][0] = ONE;  H1[1][1] = -ONE
    H2 = zero_matrix(3, 3); H2[1][1] = ONE;  H2[2][2] = -ONE
    gens.append(H1)
    gens.append(H2)
    return gens


def lift_sl3_to_gl_VW(X_sl3):
    """
    X ∈ sl(3) → X_V ∈ gl(7), the Lie-algebra lift of the SL(W)-action
    on V_W = R ⊕ W ⊕ W*.
    At the group level: g = diag(1, A, A⁻ᵀ).
    Differentiating at the identity gives X_V = diag(0, X, -Xᵀ).
    """
    X_V = zero_matrix(DIM, DIM)
    for i in range(3):
        for j in range(3):
            X_V[1 + i][1 + j] = X_sl3[i][j]
            X_V[4 + i][4 + j] = -X_sl3[j][i]
    return X_V


# ============================================================
# Action of gl(7) on Λ³(R⁷) and infinitesimal stabilizer
# ============================================================

def action_on_lambda3(X_gl7):
    """
    35×35 matrix M of the action of X ∈ gl(7) on Λ³(R⁷*),
    in the STANDARD Lie-derivative convention.

      M[row][col] = coefficient of e_{basis[row]} in X · e_{basis[col]}.

    For a 3-form ω stored as a vector w (w[col] = ω_{basis[col]}), we have
      (M · w)[row] = (X · ω)_{basis[row]}
                   = Σ_l X[l][row[0]]·ω(l, row[1], row[2]) + (cyclic).

    Implementation note. We iterate over the SOURCE basis (col), so for each
    source triple (i,j,k) and substituting position `slot` with value `l` we
    obtain a target row = sorted_index([i,j,k] with [slot] := l). The matrix
    entry needed to match the standard convention X[summed][form_position]
    is X[source[slot]][l]: the summed index is source[slot] (the value that
    left the source on its way to becoming target's missing index), and the
    matching index is l (the new value that entered).
    """
    basis = sorted_triples(DIM)
    size = len(basis)
    basis_index = {triple: idx for idx, triple in enumerate(basis)}
    M = zero_matrix(size, size)

    for col, (i, j, k) in enumerate(basis):
        contributions = [ZERO] * size
        for slot, source_index in enumerate((i, j, k)):
            for l in range(DIM):
                coef = X_gl7[source_index][l]    # standard Lie convention
                if coef == ZERO:
                    continue
                replaced = [i, j, k]
                replaced[slot] = l
                if len(set(replaced)) != 3:
                    continue          # duplicate index ⇒ vanishes
                target = basis_index[tuple(sorted(replaced))]
                contributions[target] += coef * permutation_sign(replaced)

        for row in range(size):
            M[row][col] = contributions[row]

    return M


def infinitesimal_stabilizer_constraints(omega, dim):
    """
    Build the (#triples) × (dim²) constraint matrix for X·ω = 0, X ∈ gl(dim).

    Row index r ↔ sorted triple (i,j,k).
    Col index c ↔ pair (a,b) with c = a·dim + b (matrix entry X[a][b]).
    Entry encodes the contribution of X[a][b] to (X·ω)(i,j,k).
    """
    triples = sorted_triples(dim)
    matrix = zero_matrix(len(triples), dim * dim)

    def value_at(i, j, k):
        return omega.get((i, j, k), ZERO)

    for row_idx, (i, j, k) in enumerate(triples):
        for col_idx in range(dim * dim):
            from_index = col_idx // dim    # gl-matrix row "a"
            to_index   = col_idx % dim     # gl-matrix col "b"
            entry = ZERO
            if to_index == i: entry += value_at(from_index, j, k)
            if to_index == j: entry += value_at(i, from_index, k)
            if to_index == k: entry += value_at(i, j, from_index)
            matrix[row_idx][col_idx] = entry

    return matrix


def stabilizer_dim(omega, dim):
    """
    dim of {X ∈ gl(dim) : X·ω = 0}, computed exactly.
    Returns (stabilizer_dim, rank_of_constraints).
    """
    constraints = infinitesimal_stabilizer_constraints(omega, dim)
    rank = row_reduce(constraints, dim * dim)
    return dim * dim - rank, rank


# ============================================================
# 3-form Ω_{a,b,c} on V_W, Hitchin form, and pullback
# ============================================================

def build_omega_abc(coef_pairing, coef_volW, coef_volWstar):
    """
    Ω_{a,b,c} on V_W = R⊕W⊕W*, with coordinates
      e₀ = τ,  e₁,e₂,e₃ ∈ W,  e₄,e₅,e₆ ∈ W*.

      Ω = a · τ ∧ (e₁∧e₄ + e₂∧e₅ + e₃∧e₆)   (the canonical pairing ev)
          + b · e₁∧e₂∧e₃                      (vol_W)
          + c · e₄∧e₅∧e₆                      (vol_{W*})

    Coefficients may be int / Fraction / Poly; arithmetic propagates
    transparently.
    """
    signed_sorted = {
        (0, 1, 4): coef_pairing, (0, 2, 5): coef_pairing, (0, 3, 6): coef_pairing,
        (1, 2, 3): coef_volW,
        (4, 5, 6): coef_volWstar,
    }
    return antisymmetric_form(signed_sorted)


def hitchin_form_b_over_6(omega):
    """
    Returns the matrix  b_Ω / 6,  where

      b_Ω(X, Y) · vol = (ι_X Ω) ∧ (ι_Y Ω) ∧ Ω

    is the standard Hitchin quadratic invariant (one common normalization).
    The /6 is convenient because the resulting matrix has small rational
    entries; det(b_Ω) = 6⁷ · det(b_Ω/6).

    In components:
      (b_Ω/6)[row][col] = (1/6) · Σ  sign(partition) ·
                          ω(row, γ₁, γ₂) · ω(col, γ₃, γ₄) · ω(γ₅, γ₆, γ₇),
    summed over ordered partitions of {0,…,6} into blocks
    (γ₁γ₂)(γ₃γ₄)(γ₅γ₆γ₇).
    """
    def value_at(i, j, k):
        return omega.get((i, j, k), ZERO)

    all_indices = list(range(DIM))
    b_matrix = zero_matrix(DIM, DIM)

    for row in range(DIM):
        for col in range(DIM):
            total = ZERO
            for pair_first in combinations(all_indices, 2):
                rest_after_first = [x for x in all_indices if x not in pair_first]
                for pair_second in combinations(rest_after_first, 2):
                    last_triple = tuple(x for x in rest_after_first
                                        if x not in pair_second)
                    if len(last_triple) != 3:
                        continue

                    val_row    = value_at(row, pair_first[0], pair_first[1])
                    val_col    = value_at(col, pair_second[0], pair_second[1])
                    val_triple = value_at(last_triple[0], last_triple[1], last_triple[2])
                    if val_row == ZERO or val_col == ZERO or val_triple == ZERO:
                        continue

                    full_perm = list(pair_first) + list(pair_second) + list(last_triple)
                    if len(set(full_perm)) != DIM:
                        continue

                    total += permutation_sign(full_perm) * val_row * val_col * val_triple

            b_matrix[row][col] = total / Fraction(6)

    return b_matrix


def pullback_3form(H, omega, dim):
    """(H*ω)_{ijk} = Σ_{abc} H[a][i] H[b][j] H[c][k] · ω_{abc}."""
    pulled = {}
    for i in range(dim):
        for j in range(dim):
            for k in range(dim):
                val = ZERO
                for a in range(dim):
                    for b in range(dim):
                        for c in range(dim):
                            coef = omega.get((a, b, c), ZERO)
                            if coef == ZERO:
                                continue
                            val += H[a][i] * H[b][j] * H[c][k] * coef
                pulled[(i, j, k)] = val
    return pulled


# ============================================================
# PART 1: dim SL(W)-invariant 3-forms = 3
# ============================================================

print("=" * 60)
print("PART 1: SL(W)-invariant 3-forms on V_W")
print("=" * 60)

# Stack the action matrices for all sl(3) generators.
# ω is SL(W)-invariant iff (stacked matrix) · ω = 0.
stacked_constraints = []
for X_sl3 in sl3_generators():
    X_V = lift_sl3_to_gl_VW(X_sl3)
    action_matrix = action_on_lambda3(X_V)
    for row in action_matrix:
        stacked_constraints.append(row)

invariant_dim = kernel_dimension(stacked_constraints, 35)
print(f"  dim SL(W)-invariant 3-forms = {invariant_dim}")
require(invariant_dim == 3, f"Expected 3, got {invariant_dim}")
print("  ✓ VERIFIED: exactly 3 invariants.\n")


# ============================================================
# PART 2: SYMBOLIC verification of b_Ω/6 block structure
# ============================================================

print("=" * 60)
print("PART 2: Symbolic proof of det(b_Ω) = 4374·a⁹b⁶c⁶")
print("=" * 60)

# Build Ω with symbolic a, b, c ∈ Q[a,b,c]; run hitchin_form_b_over_6 on it;
# check every entry of the resulting matrix against the predicted formula
# AS POLYNOMIALS. Equality of polynomials in Q[a,b,c] is a definitive
# identity (no sampling, no Schwartz-Zippel — direct comparison of
# coefficient dicts).
#
# Predicted (b_Ω/6) block structure:
#   (b/6)[0,0]      = -a³
#   (b/6)[i,i+3]    = (b/6)[i+3,i] = abc/2     for i = 1, 2, 3
#   all other entries = 0

a_sym = Poly.variable(0)
b_sym = Poly.variable(1)
c_sym = Poly.variable(2)

Omega_sym = build_omega_abc(a_sym, b_sym, c_sym)
b_over_6_sym = hitchin_form_b_over_6(Omega_sym)

expected_diag   = -(a_sym ** 3)
expected_offdiag = (a_sym * b_sym * c_sym) / 2

structure_ok = True
mismatches = []

for r in range(DIM):
    for c in range(DIM):
        if r == 0 and c == 0:
            target = expected_diag
        elif 1 <= r <= 3 and c == r + 3:
            target = expected_offdiag
        elif 4 <= r <= 6 and c == r - 3:
            target = expected_offdiag
        else:
            target = Poly()    # zero polynomial

        if b_over_6_sym[r][c] != target:
            structure_ok = False
            mismatches.append((r, c, b_over_6_sym[r][c], target))

if structure_ok:
    print("  Symbolic check: every entry of b_Ω/6 in Q[a,b,c] matches the")
    print(f"  predicted block structure (49 polynomial identities, all hold).")
    print(f"  (b/6)[0][0]          = {b_over_6_sym[0][0]!r}        (expected: {expected_diag!r})")
    print(f"  (b/6)[1][4]          = {b_over_6_sym[1][4]!r}    (expected: {expected_offdiag!r})")
    print(f"  (b/6)[2][5]          = {b_over_6_sym[2][5]!r}    (expected: {expected_offdiag!r})")
    print(f"  (b/6)[3][6]          = {b_over_6_sym[3][6]!r}    (expected: {expected_offdiag!r})")
    print(f"  ✓ VERIFIED: block structure proven symbolically over Q[a,b,c].")
else:
    for (r, c, got, want) in mismatches[:5]:
        print(f"  ✗ Mismatch at [{r}][{c}]: got {got!r}, expected {want!r}")
    require(False, f"Symbolic block-structure verification failed "
                   f"({len(mismatches)} mismatches).")

# Direct symbolic determinant check. This is redundant once the block structure
# has been proven, but it makes the determinant identity explicit in the
# verifier output.
det_b6_sym_direct = exact_det(b_over_6_sym, DIM)
expected_det_b6_sym = (a_sym ** 9) * (b_sym ** 6) * (c_sym ** 6) / 64
require(det_b6_sym_direct == expected_det_b6_sym,
        f"Expected det(b_Ω/6)={expected_det_b6_sym!r}, "
        f"got {det_b6_sym_direct!r}")
print(f"  Direct symbolic determinant check: det(b_Ω/6) = {det_b6_sym_direct!r} ✓")

# Det formula also follows transparently from the block structure:
#   det(b/6) = (-a³) · ((abc/2) · (abc/2))³ · sign of off-diag block permutation
#            = (-a³) · (a²b²c²/4)³ · (-1)³           [each anti-diag pair → −]
# Carefully:
#   The matrix splits as 1×1 block (−a³) and three 2×2 blocks of the form
#     [[0, abc/2], [abc/2, 0]],  each with det = −a²b²c²/4.
#   Total: (−a³) · (−a²b²c²/4)³ = a⁹b⁶c⁶ / 64.
# Hence det(b_Ω) = 6⁷ · a⁹b⁶c⁶ / 64 = 4374 · a⁹b⁶c⁶.
print()
print("  Consequence (block-structure determinant of a 1+2+2+2 block matrix):")
print("     det(b_Ω/6) = (-a³) · (-a²b²c²/4)³ = a⁹b⁶c⁶/64")
print("     det(b_Ω)   = 6⁷ · a⁹b⁶c⁶/64       = 4374·a⁹b⁶c⁶")
print()

# Numerical sanity demonstrations (illustrative, not part of the proof):
print("  Numerical illustrations:")
for (alpha, beta, gamma) in [(1, 1, 1), (2, 3, 7), (-2, 3, -5)]:
    Omega_num = build_omega_abc(alpha, beta, gamma)
    b_num = hitchin_form_b_over_6(Omega_num)
    det_b6 = Fraction(alpha**9 * beta**6 * gamma**6, 64)
    det_b  = det_b6 * 6**7
    require(b_num[0][0] == Fraction(-alpha**3), "block diag mismatch")
    require(b_num[1][4] == Fraction(alpha*beta*gamma, 2), "block offdiag mismatch")
    print(f"    (a,b,c)=({alpha},{beta},{gamma}): "
          f"det(b_Ω/6) = {det_b6}, det(b_Ω) = {det_b}  ✓")
print()


# ============================================================
# PART 3: Ω_{1,1,1} ⇒ sig(b_Ω) = (3,4) and dim Stab = 14
# ============================================================

print("=" * 60)
print("PART 3: Ω_{1,1,1} — sig(b_Ω) = (3,4), dim Stab = 14")
print("=" * 60)

Omega_111 = build_omega_abc(1, 1, 1)
b_111     = hitchin_form_b_over_6(Omega_111)

# Sanity from the block structure (a = b = c = 1):
#   b/6 = diag(-1) ⊕ [[0, 1/2], [1/2, 0]]³
#   Each off-diagonal block has eigenvalues ±1/2 → signature (1,1).
#   Three blocks → (3,3). Plus the −1 entry → total (3,4).
# Multiplying by 6 > 0 preserves signature, so sig(b_Ω) = sig(b_Ω/6) = (3,4).
pos, neg = exact_signature(b_111, DIM)
print(f"  sig(b_Ω) = sig(b_Ω/6) = ({pos},{neg}) "
      f"(exact congruence diagonalization)")
require((pos, neg) == (3, 4), f"Expected (3,4), got ({pos},{neg})")
print("  ✓ sig = (3,4). Exact.")
print("  Note: dim Stab = 14 and sig = (3,4) are the numerical invariants of")
print("        G₂_split. Identifying the stabilizer with G₂_split itself requires")
print("        the standard classification theorem for stable 3-forms in dimension 7;")
print("        this script verifies the computational invariants, not that external theorem.\n")

stab_dim_VW, _ = stabilizer_dim(Omega_111, DIM)
print(f"  dim Stab(Ω_{{1,1,1}}) = {stab_dim_VW} (exact)")
require(stab_dim_VW == 14, f"Expected 14, got {stab_dim_VW}")
print()


# ============================================================
# PART 4: Ω_W ~ Ω_old as 3-forms on R⁷
# ============================================================

print("=" * 60)
print("PART 4: Ω_W ~ Ω_old as 3-forms on R⁷")
print("=" * 60)

# We exhibit an integer matrix H with det(H) = 8 (so H ∈ GL(7,Q), NOT GL(7,Z))
# satisfying  H*Ω_W = 2·Ω_old.
# To get H'*Ω_W = Ω_old exactly over R, rescale: H' := 2^(-1/3) · H. Then
#   H'*Ω_W = (2^(-1/3))³ · H*Ω_W = (1/2) · 2 · Ω_old = Ω_old.
#
# The fact that this displayed rescaling is irrational is not by itself a proof
# of non-equivalence over Q: a different rational witness might have existed.
# We therefore add a determinant obstruction below. For any K ∈ GL(7,Q),
# covariance of the Hitchin bilinear form gives
#   B_{K*Ω} = det(K) · Kᵀ B_Ω K,
# hence
#   det(B_{K*Ω}) = det(K)^9 det(B_Ω).
# Thus exact rational equivalence K*Ω_W = Ω_old would force
#   det(B_old) / det(B_W) = det(K)^9,
# i.e. the determinant ratio must be a ninth power in Q.

H_witness = [
    [Fraction(0), Fraction(0), Fraction(0), Fraction( 1), Fraction(0), Fraction(0), Fraction( 0)],
    [Fraction(1), Fraction(0), Fraction(0), Fraction( 0), Fraction(-1), Fraction(0), Fraction( 0)],
    [Fraction(0), Fraction(1), Fraction(0), Fraction( 0), Fraction(0), Fraction(-1), Fraction( 0)],
    [Fraction(0), Fraction(0), Fraction(1), Fraction( 0), Fraction(0), Fraction(0), Fraction( 1)],
    [Fraction(1), Fraction(0), Fraction(0), Fraction( 0), Fraction(1), Fraction(0), Fraction( 0)],
    [Fraction(0), Fraction(1), Fraction(0), Fraction( 0), Fraction(0), Fraction(1), Fraction( 0)],
    [Fraction(0), Fraction(0), Fraction(1), Fraction( 0), Fraction(0), Fraction(0), Fraction(-1)],
]

Omega_W = antisymmetric_form({
    (0, 1, 4): +1, (0, 2, 5): +1, (0, 3, 6): +1,
    (1, 2, 3): +1,
    (4, 5, 6): +1,
})
# Ω_old: classical 7-term Fano-plane representative (used in proofs 22–27).
# In the same GL(7,R)-orbit as Ω_W.
Omega_old = antisymmetric_form({
    (0, 1, 2): +1, (0, 3, 4): -1, (0, 5, 6): -1,
    (1, 3, 5): -1, (2, 4, 5): +1,
    (1, 4, 6): +1, (2, 3, 6): +1,
})

pulled = pullback_3form(H_witness, Omega_W, DIM)
max_diff = ZERO
for i in range(DIM):
    for j in range(DIM):
        for k in range(DIM):
            expected = Omega_old.get((i, j, k), ZERO) * 2
            actual   = pulled.get((i, j, k),    ZERO)
            diff = abs(actual - expected)
            if diff > max_diff:
                max_diff = diff

print(f"  Exact witness H (integer entries): max|H*Ω_W − 2·Ω_old| = {max_diff}")

det_H = exact_det(H_witness, DIM)
print(f"  det(H) = {det_H}  →  H ∈ GL(7,Q) ⊂ GL(7,R),  H ∉ GL(7,Z)")
require(det_H != ZERO, "det(H) = 0, not invertible!")
require(max_diff == ZERO, f"Exact witness failed! diff={max_diff}")

B_W_over_6 = hitchin_form_b_over_6(Omega_W)
B_old_over_6 = hitchin_form_b_over_6(Omega_old)
det_B_W_over_6 = exact_det(B_W_over_6, DIM)
det_B_old_over_6 = exact_det(B_old_over_6, DIM)
det_ratio = det_B_old_over_6 / det_B_W_over_6

print(f"  det(B_ΩW/6)   = {det_B_W_over_6}")
print(f"  det(B_Ωold/6) = {det_B_old_over_6}")
print(f"  determinant ratio det(B_old)/det(B_W) = {det_ratio}")
require(det_ratio == Fraction(64), f"Expected determinant ratio 64, got {det_ratio}")
require(not is_rational_nth_power(det_ratio, 9),
        "Unexpected: determinant ratio 64 is a ninth power in Q")
print("  Determinant obstruction over Q: 64 is not a ninth power in Q.")
print("  Equivalently: v_2(q^9) ∈ 9Z for rational q, but v_2(64) = 6 ∉ 9Z.")
print("  Therefore no K ∈ GL(7,Q) can satisfy K*Ω_W = Ω_old exactly.")
print("  Rescaling H' = 2^(-1/3)·H ∈ GL(7,R) gives H'*Ω_W = Ω_old exactly.")
print("  ✓ VERIFIED: Ω_W and Ω_old lie in the same GL(7,R)-orbit, not the same")
print("    exact GL(7,Q)-orbit. Over Q the witness gives Ω_W ~ 2·Ω_old.\n")


# ============================================================
# PART 5: dim=6 (W⊕W*) → classical stabilizer (not G₂)
# ============================================================

print("=" * 60)
print("PART 5: On W⊕W* (dim 6), stabilizer is classical")
print("=" * 60)

# Ω₆ = vol_W + vol_{W*} on R⁶, coordinates 0,1,2 = W and 3,4,5 = W*.
Omega_6 = antisymmetric_form({
    (0, 1, 2): +1,
    (3, 4, 5): +1,
})

stab_dim_6, rank_6 = stabilizer_dim(Omega_6, 6)
print(f"  Ω₆ = vol_W + vol_{{W*}} on R⁶")
print(f"  rank = {rank_6}, dim Stab = {stab_dim_6}")
require(stab_dim_6 == 16, f"Expected 16, got {stab_dim_6}")
print(f"  ✓ dim Stab = {stab_dim_6} > 14 → NOT G₂ (classical).\n")


# ============================================================
# PART 6: SL(W) preserves Ω_W (exact algebraic proof)
# ============================================================

print("=" * 60)
print("PART 6: SL(W) preserves Ω_W")
print("=" * 60)

# Exact proof. For A ∈ SL(W), the action on V_W = R⊕W⊕W* is
#   g = diag(1, A, A⁻ᵀ).
# Then:
#   g*(τ ∧ ev)    = τ ∧ ev
#                  ← The pairing ev ∈ W* ⊗ W is invariant because the
#                    actions on W (by A) and on W* (by A⁻ᵀ) are
#                    contragredient: (A⁻ᵀφ)(Aw) = φ(A⁻¹·A·w) = φ(w).
#   g*(vol_W)     = det(A) · vol_W      = vol_W      (det A = 1)
#   g*(vol_{W*})  = det(A⁻ᵀ) · vol_{W*} = vol_{W*}   (det A⁻ᵀ = 1)
# Hence g*(Ω_W) = Ω_W. QED.
#
# Computational check on all 8 sl(3) generators (infinitesimal form):
# X · Ω_W must vanish identically.

all_zero = True
generators_count = 0
for gen_idx, X_sl3 in enumerate(sl3_generators()):
    generators_count += 1
    X_V = lift_sl3_to_gl_VW(X_sl3)

    for (i, j, k) in sorted_triples(DIM):
        val = ZERO
        for l in range(DIM):
            val += (X_V[l][i] * Omega_W.get((l, j, k), ZERO)
                  + X_V[l][j] * Omega_W.get((i, l, k), ZERO)
                  + X_V[l][k] * Omega_W.get((i, j, l), ZERO))
        if val != ZERO:
            all_zero = False
            print(f"  ✗ Generator {gen_idx}: X·Ω[{i},{j},{k}] = {val}")
            break

require(all_zero, "SL(W) does not preserve Ω_W!")
print(f"  All {generators_count} sl(3) generators annihilate Ω_W.")
print("  ✓ VERIFIED: SL(W) preserves Ω_W. Exact (Fraction).\n")


# ============================================================
# SUMMARY
# ============================================================

print("=" * 60)
print("SUMMARY")
print("=" * 60)
print("  1. dim SL(W)-invariant 3-forms = 3                       ✓")
print("  2. det(b_Ω) = 4374·a⁹b⁶c⁶  (symbolic proof in Q[a,b,c])  ✓")
print("  3. Ω_{1,1,1}: dim Stab = 14, sig(b_Ω) = (3,4)            ✓")
print("     (G₂_split identification uses a standard external classification theorem)")
print("  4. Ω_W ~ Ω_old over R; NOT GL(7,Q)-equivalent, since det-ratio 64 ∉ Q^9 ✓")
print("  5. dim=6: dim Stab = 16 > 14 → not G₂ (classical)        ✓")
print("  6. SL(W) preserves Ω_W (all 8 sl(3) generators)          ✓")
print("  All internal computational claims of Proof 21 verified.")
print("  External G₂_split identification rests on the standard classification/stable-form theorem.")
tweet_url

    
SHA-256
Visualization