#!/usr/bin/env python3
"""Verify the five-simple quiver algebra using exact path-space arithmetic.

The input is the ten-arrow quiver and the quadratic relations displayed in
main.tex (Theorem labelled thm:c3-quiver).  This verifier is independent of the link-state/spin-chain
matrix reconstruction.  It uses only the Python standard library and performs
Gaussian elimination over the rational numbers, with no floating-point or
finite-field calculation.

Paths are tuples of vertices in traversal order.  Thus ('a', 'b', 'd')
represents u_bd u_ab in the manuscript's rightmost-arrow-first convention.
Cartan rows index simple modules and columns index projective covers.  Ext^2
rows index the first argument and columns index the second argument.
"""

from fractions import Fraction


VERTICES = tuple("abcde")
ARROWS = (
    ("a", "b"), ("a", "c"),
    ("b", "a"), ("b", "d"), ("b", "e"),
    ("c", "d"), ("c", "e"),
    ("d", "b"), ("d", "c"),
    ("e", "c"),
)

# Within each group all displayed paths are equated.  Every other quadratic
# path is zero.  The singleton group is the one additional nonzero path.
NONZERO_QUADRATIC_GROUPS = (
    (("a", "b", "d"), ("a", "c", "d")),
    (("a", "b", "e"), ("a", "c", "e")),
    (("b", "a", "c"), ("b", "d", "c"), ("b", "e", "c")),
    (("d", "b", "e"),),
)

PHYSICAL_SIMPLE_DIMENSIONS = (1, 13, 4, 1, 1)

EXPECTED_CARTAN = (
    (1, 1, 0, 0, 0),
    (1, 1, 0, 1, 0),
    (1, 1, 1, 1, 1),
    (1, 1, 1, 1, 0),
    (1, 1, 1, 1, 1),
)

EXPECTED_EXT2 = (
    (1, 0, 0, 1, 1),
    (0, 2, 2, 0, 0),
    (0, 1, 2, 0, 0),
    (1, 0, 0, 2, 1),
    (0, 0, 0, 1, 1),
)

EXPECTED_PROJECTIVE_DIMENSIONS = (20, 20, 6, 19, 5)


def require(condition, message):
    """Do not disable mathematical checks when Python is run with -O."""
    if not condition:
        raise RuntimeError("VERIFICATION FAILED: " + message)


def enumerate_paths(length):
    paths = [(vertex,) for vertex in VERTICES]
    for _ in range(length):
        paths = [path + (target,)
                 for path in paths
                 for source, target in ARROWS
                 if source == path[-1]]
    return tuple(sorted(paths))


def rank(rows, width):
    """Exact rational rank, including the empty-matrix cases."""
    if width == 0:
        return 0
    echelon = {}
    for row in rows:
        require(len(row) == width, "inconsistent matrix row width")
        vector = [Fraction(entry) for entry in row]
        for pivot in sorted(echelon):
            coefficient = vector[pivot]
            if coefficient:
                basis_vector = echelon[pivot]
                vector = [x - coefficient*y
                          for x, y in zip(vector, basis_vector)]
        pivot = next((i for i, value in enumerate(vector) if value), None)
        if pivot is not None:
            coefficient = vector[pivot]
            echelon[pivot] = [value/coefficient for value in vector]
    return len(echelon)


def dense(relation, basis):
    return [relation.get(path, 0) for path in basis]


def quadratic_relations(paths2):
    nonzero = {path for group in NONZERO_QUADRATIC_GROUPS for path in group}
    require(nonzero <= set(paths2), "a displayed path is not composable")
    relations = [{path: 1} for path in paths2 if path not in nonzero]
    zero_relations = len(relations)
    for group in NONZERO_QUADRATIC_GROUPS:
        first = group[0]
        for path in group[1:]:
            require((first[0], first[-1]) == (path[0], path[-1]),
                    "a quadratic relation is not homogeneous in endpoints")
            relations.append({path: 1, first: -1})
    return relations, zero_relations


def degree_three_ideal(relations2):
    """Generate I_3 = Q_1 I_2 + I_2 Q_1, with every compatible arrow."""
    relations3 = []
    for relation in relations2:
        endpoints = {(path[0], path[-1]) for path in relation}
        require(len(endpoints) == 1, "inhomogeneous relation endpoints")
        source, target = next(iter(endpoints))
        for tail, head in ARROWS:
            if head == source:
                relations3.append({(tail,) + path: coefficient
                                   for path, coefficient in relation.items()})
            if tail == target:
                relations3.append({path + (head,): coefficient
                                   for path, coefficient in relation.items()})
    return relations3


def print_matrix(title, matrix):
    print(title)
    for label, row in zip(VERTICES, matrix):
        print("  " + label + ": " + " ".join(str(entry) for entry in row))


def main():
    paths0 = enumerate_paths(0)
    paths1 = enumerate_paths(1)
    paths2 = enumerate_paths(2)
    paths3 = enumerate_paths(3)
    require(len(paths1) == 10, "the quiver must have ten arrows")
    require(len(paths2) == 20, "the quadratic path count must be twenty")
    require(len(paths3) == 38, "the cubic path count must be thirty-eight")

    relations2, zero_count = quadratic_relations(paths2)
    rows2 = [dense(relation, paths2) for relation in relations2]
    rank2 = rank(rows2, len(paths2))
    require(zero_count == 12, "there must be twelve zero-path relations")
    require(len(relations2) == 16 and rank2 == 16,
            "the quadratic relations must span a sixteen-dimensional space")
    require(len(paths2) - rank2 == 4,
            "the degree-two quotient must be four-dimensional")

    relations3 = degree_three_ideal(relations2)
    require(all(set(relation) <= set(paths3) for relation in relations3),
            "an ideal generator contains a noncomposable cubic path")
    rank3 = rank([dense(relation, paths3) for relation in relations3], len(paths3))
    require(rank3 == len(paths3),
            "the quadratic relations must annihilate every cubic path")

    # Resolve every homogeneous endpoint corner separately.  Since the ideal
    # has no terms of degree zero or one, its degree-two components are the
    # minimal relation spaces computing Ext^2 between the corresponding simples.
    second_layer = [[0 for _ in VERTICES] for _ in VERTICES]
    ext2 = [[0 for _ in VERTICES] for _ in VERTICES]
    cartan = [[0 for _ in VERTICES] for _ in VERTICES]
    first_layer = [[0 for _ in VERTICES] for _ in VERTICES]
    for source_index, source in enumerate(VERTICES):
        for target_index, target in enumerate(VERTICES):
            corner = tuple(path for path in paths2
                           if path[0] == source and path[-1] == target)
            relation_rank = rank([dense(relation, corner)
                                  for relation in relations2], len(corner))
            quotient_dimension = len(corner) - relation_rank
            arrow_count = sum(path == (source, target) for path in paths1)
            second_layer[target_index][source_index] = quotient_dimension
            first_layer[target_index][source_index] = arrow_count
            ext2[source_index][target_index] = relation_rank
            cartan[target_index][source_index] = (
                int(source == target) + arrow_count + quotient_dimension)

    cartan = tuple(tuple(row) for row in cartan)
    ext2 = tuple(tuple(row) for row in ext2)
    require(cartan == EXPECTED_CARTAN, "Cartan matrix differs from the manuscript")
    require(ext2 == EXPECTED_EXT2, "Ext^2 matrix differs from the manuscript")
    require(sum(map(sum, cartan)) == 19, "the basic algebra dimension must be 19")
    require(sum(map(sum, ext2)) == 16, "the total minimal-relation count must be 16")

    sizes = PHYSICAL_SIMPLE_DIMENSIONS
    projective_dimensions = tuple(
        sum(sizes[i]*cartan[i][j] for i in range(len(VERTICES)))
        for j in range(len(VERTICES)))
    require(projective_dimensions == EXPECTED_PROJECTIVE_DIMENSIONS,
            "physical projective dimensions differ from the manuscript")
    semisimple_dimension = sum(value*value for value in sizes)
    radical_layer_dimension = sum(
        sizes[i]*first_layer[i][j]*sizes[j]
        for i in range(len(VERTICES)) for j in range(len(VERTICES)))
    radical_square_dimension = sum(
        sizes[i]*second_layer[i][j]*sizes[j]
        for i in range(len(VERTICES)) for j in range(len(VERTICES)))
    full_dimension = sum(sizes[j]*projective_dimensions[j]
                         for j in range(len(VERTICES)))
    require((semisimple_dimension, radical_layer_dimension,
             radical_square_dimension, full_dimension) == (188, 85, 55, 328),
            "physical layer dimensions differ from the manuscript")

    print("Exact quiver verification over the rational numbers")
    print("Vertices: a b c d e; arrows: 10")
    print(f"Degree-two paths: {len(paths2)}; relation rank: {rank2}; "
          f"quotient dimension: {len(paths2)-rank2}")
    print(f"Degree-three paths: {len(paths3)}; generated ideal rank: {rank3}; "
          f"quotient dimension: {len(paths3)-rank3}")
    print(f"Cubic ideal generators checked: {len(relations3)}")
    print("Basic algebra dimension: 5 + 10 + 4 = 19; all higher degrees vanish.")
    print_matrix("Cartan matrix (rows: simples; columns: projectives):", cartan)
    print("Physical projective dimensions:", " ".join(map(str, projective_dimensions)))
    print("Physical algebra layers: 188 + 85 + 55 = 328")
    print_matrix("Ext^2 matrix (rows: first argument; columns: second argument):", ext2)
    print("Cartan rank:", rank(cartan, len(VERTICES)))
    print("PASS: path spaces, cubic vanishing, Cartan matrix, projective dimensions, "
          "and Ext^2 matrix agree with the manuscript.")


if __name__ == "__main__":
    main()