verify-secant-gap-boundary.py

#!/usr/bin/env python3
"""Exact audit for the Sylvester inverse-system rank gap.

The script verifies two claims over Q.

1. A concrete three-plane W of zero-diagonal quadrics in four variables
   contains no nonzero quadric of rank at most two.
2. If I is generated by W^perp in degree two, then S/I has
   beta_(2,4)=3, so I is not linearly presented.

It also runs the existing exact flat-Koszul engine on the full PG(3,2)
positive control.  No external algebra package is required.
"""

from __future__ import annotations

import importlib.util
import sys
from fractions import Fraction
from itertools import combinations
from pathlib import Path


EDGES = list(combinations(range(4), 2))
EDGE_INDEX = {edge: index for index, edge in enumerate(EDGES)}

# Rows are indexed by 12,13,14,23,24,34.  A parameter vector u in Q^3
# gives the six off-diagonal coefficients C u of a quadric in W.
C = [
    [-1, -1, -1],
    [1, -1, -1],
    [0, 1, -1],
    [0, 0, -1],
    [-1, 0, 1],
    [0, 1, 1],
]


def determinant_3(matrix: list[list[int]]) -> int:
    return (
        matrix[0][0]
        * (matrix[1][1] * matrix[2][2] - matrix[1][2] * matrix[2][1])
        - matrix[0][1]
        * (matrix[1][0] * matrix[2][2] - matrix[1][2] * matrix[2][0])
        + matrix[0][2]
        * (matrix[1][0] * matrix[2][1] - matrix[1][1] * matrix[2][0])
    )


def cross(left: list[int], right: list[int]) -> list[int]:
    return [
        left[1] * right[2] - left[2] * right[1],
        left[2] * right[0] - left[0] * right[2],
        left[0] * right[1] - left[1] * right[0],
    ]


def dot(left: list[int], right: list[int]) -> int:
    return sum(a * b for a, b in zip(left, right))


def rref(matrix: list[list[int | Fraction]]) -> tuple[list[list[Fraction]], list[int]]:
    work = [[Fraction(entry) for entry in row] for row in matrix]
    row_count = len(work)
    column_count = len(work[0]) if work else 0
    pivot_columns: list[int] = []
    pivot_row = 0
    for column in range(column_count):
        pivot = next(
            (row for row in range(pivot_row, row_count) if work[row][column]),
            None,
        )
        if pivot is None:
            continue
        work[pivot_row], work[pivot] = work[pivot], work[pivot_row]
        leading = work[pivot_row][column]
        work[pivot_row] = [entry / leading for entry in work[pivot_row]]
        for row in range(row_count):
            if row == pivot_row or not work[row][column]:
                continue
            factor = work[row][column]
            work[row] = [
                entry - factor * pivot_entry
                for entry, pivot_entry in zip(work[row], work[pivot_row])
            ]
        pivot_columns.append(column)
        pivot_row += 1
        if pivot_row == row_count:
            break
    return work, pivot_columns


def rank(matrix: list[list[int | Fraction]]) -> int:
    return len(rref(matrix)[1])


def nullspace(matrix: list[list[int]]) -> list[list[Fraction]]:
    reduced, pivots = rref(matrix)
    column_count = len(matrix[0])
    free_columns = [column for column in range(column_count) if column not in pivots]
    basis: list[list[Fraction]] = []
    for free in free_columns:
        vector = [Fraction(0)] * column_count
        vector[free] = Fraction(1)
        for row, pivot in enumerate(pivots):
            vector[pivot] = -reduced[row][free]
        basis.append(vector)
    return basis


def audit_rank_gap() -> None:
    singleton_determinants = []
    for singleton in range(4):
        triple = [element for element in range(4) if element != singleton]
        rows = [C[EDGE_INDEX[edge]] for edge in combinations(triple, 2)]
        singleton_determinants.append(determinant_3(rows))
    assert singleton_determinants == [1, 2, -3, -2]

    balanced_determinants = []
    kernel_vectors = []
    for left in [(0, 1), (0, 2), (0, 3)]:
        right = tuple(element for element in range(4) if element not in left)
        kernel = cross(
            C[EDGE_INDEX[tuple(sorted(left))]],
            C[EDGE_INDEX[tuple(sorted(right))]],
        )
        cross_block = [
            [
                dot(C[EDGE_INDEX[tuple(sorted((a, b)))]], kernel)
                for b in right
            ]
            for a in left
        ]
        determinant = (
            cross_block[0][0] * cross_block[1][1]
            - cross_block[0][1] * cross_block[1][0]
        )
        kernel_vectors.append(kernel)
        balanced_determinants.append(determinant)
    assert kernel_vectors == [[0, 1, -1], [-1, 0, -1], [-1, 0, 0]]
    assert balanced_determinants == [-2, -3, 1]

    print("rank-gap cut certificate")
    print("  singleton cut determinants:", singleton_determinants)
    print("  balanced-cut kernel vectors:", kernel_vectors)
    print("  balanced-cut cross determinants:", balanced_determinants)
    print("  conclusion: P(W) contains no rank-one or rank-two quadric")


def audit_nonlinear_syzygies() -> None:
    # The transpose kernel gives the three nonsquare generators of W^perp.
    relation_basis = nullspace(
        [[C[edge][coordinate] for edge in range(6)] for coordinate in range(3)]
    )
    expected = [
        [Fraction(-1, 4), Fraction(-1, 4), Fraction(-1, 2), 1, 0, 0],
        [Fraction(-1, 4), Fraction(3, 4), Fraction(1, 2), 0, 1, 0],
        [Fraction(1, 2), Fraction(1, 2), 0, 0, 0, 1],
    ]
    assert relation_basis == expected

    triples = list(combinations(range(4), 3))
    triple_index = {triple: index for index, triple in enumerate(triples)}
    cubic_relations: list[list[Fraction]] = []
    for relation in relation_basis:
        for variable in range(4):
            row = [Fraction(0)] * 4
            for edge, coefficient in zip(EDGES, relation):
                if variable not in edge:
                    triple = tuple(sorted(edge + (variable,)))
                    row[triple_index[triple]] += coefficient
            cubic_relations.append(row)
    assert rank(cubic_relations) == 4

    # Since A_3=0, every element of A_2 tensor wedge^2 is a cycle.
    # Compute d_3: A_1 tensor wedge^3 -> A_2 tensor wedge^2 exactly.
    wedges_2 = list(combinations(range(4), 2))
    wedge_2_index = {wedge: index for index, wedge in enumerate(wedges_2)}
    wedges_3 = list(combinations(range(4), 3))
    columns: list[list[Fraction]] = []
    for coefficient_variable in range(4):
        for wedge in wedges_3:
            column = [Fraction(0)] * 18
            for position, exterior_variable in enumerate(wedge):
                if exterior_variable == coefficient_variable:
                    continue
                pair = tuple(sorted((coefficient_variable, exterior_variable)))
                output_wedge = wedge[:position] + wedge[position + 1 :]
                wedge_index = wedge_2_index[output_wedge]
                sign = -1 if position % 2 else 1
                for coordinate, value in enumerate(C[EDGE_INDEX[pair]]):
                    column[coordinate * 6 + wedge_index] += sign * value
            columns.append(column)
    differential = [
        [columns[column][row] for column in range(16)] for row in range(18)
    ]
    differential_rank = rank(differential)
    assert differential_rank == 15
    h_2_4 = 18 - differential_rank
    assert h_2_4 == 3

    print("apolar algebra certificate")
    print("  dim A_2: 3")
    print("  rank of cubic relation span: 4 (hence A_3=0)")
    print("  dim K_(2,4): 18")
    print("  rank d_3: 15")
    print("  dim H_(2,4) = beta_(2,4):", h_2_4)


def audit_pg32_positive_control() -> None:
    here = Path(__file__).resolve().parent
    engine_path = here / "verify-flat-koszul-blocks.py"
    spec = importlib.util.spec_from_file_location("flat_koszul_blocks", engine_path)
    assert spec is not None and spec.loader is not None
    module = importlib.util.module_from_spec(spec)
    sys.modules[spec.name] = module
    spec.loader.exec_module(module)
    columns = [
        tuple((value >> bit) & 1 for bit in range(4)) for value in range(1, 16)
    ]
    matroid = module.vector_matroid("PG(3,2)", columns, 2)
    lines = module.audit(matroid)
    for degree in (4, 5, 6):
        line = next(item for item in lines if f"degree {degree}:" in item)
        assert "H2=0" in line
    print("full PG(3,2) positive control")
    print("\n".join(lines))


if __name__ == "__main__":
    audit_rank_gap()
    print()
    audit_nonlinear_syzygies()
    print()
    audit_pg32_positive_control()