Abelian categories

"""The snake lemma, computed: kernels, cokernels and the connecting map mod p.

Vector spaces over F_p form an abelian category in which every short exact
sequence splits, so every morphism of short exact sequences
0 -> A -> B -> C -> 0  to  0 -> A' -> B' -> C' -> 0 is, up to isomorphism,
B = A + C, B' = A' + C' with the middle map in block form
    beta = [[alpha, x], [0, gamma]],   x : C -> A' arbitrary.
For random such diagrams we compute ker and coker of alpha, beta, gamma,
the connecting map delta(w) = [x w] in coker(alpha), and check that the
six-term sequence
    0 -> ker a -> ker b -> ker c --delta--> coker a -> coker b -> coker c -> 0
is exact at every term (image dimension = kernel dimension, composites 0).
Exact modular arithmetic; seeded random examples.
"""

import random

P = 5

def rank(cols, p=P):
    """Rank over F_p of the matrix whose columns are the given vectors."""
    rows = [list(r) for r in zip(*cols)] if cols else []
    r = 0
    ncols = len(cols)
    for c in range(ncols):
        pivot = next((i for i in range(r, len(rows)) if rows[i][c] % p), None)
        if pivot is None:
            continue
        rows[r], rows[pivot] = rows[pivot], rows[r]
        inv = pow(rows[r][c], p - 2, p)
        rows[r] = [v * inv % p for v in rows[r]]
        for i in range(len(rows)):
            if i != r and rows[i][c] % p:
                k = rows[i][c]
                rows[i] = [(a - k * b) % p for a, b in zip(rows[i], rows[r])]
        r += 1
    return r

def apply(mat, vec, p=P):
    """Matrix (list of rows) times column vector, mod p."""
    return [sum(a * b for a, b in zip(row, vec)) % p for row in mat]

def columns(mat):
    """The columns of a matrix given as a list of rows."""
    return [list(c) for c in zip(*mat)] if mat and mat[0] else []

def identity(n):
    """The n x n identity matrix."""
    return [[int(i == j) for j in range(n)] for i in range(n)]

def kernel_basis(mat, n, p=P):
    """A basis of the kernel of mat : F_p^n -> F_p^m (mat has n columns)."""
    rows = [list(r) for r in mat]
    pivots = []
    r = 0
    for c in range(n):
        pivot = next((i for i in range(r, len(rows)) if rows[i][c] % p), None)
        if pivot is None:
            continue
        rows[r], rows[pivot] = rows[pivot], rows[r]
        inv = pow(rows[r][c], p - 2, p)
        rows[r] = [v * inv % p for v in rows[r]]
        for i in range(len(rows)):
            if i != r and rows[i][c] % p:
                k = rows[i][c]
                rows[i] = [(a - k * b) % p for a, b in zip(rows[i], rows[r])]
        pivots.append(c)
        r += 1
    basis = []
    for free in (c for c in range(n) if c not in pivots):
        v = [0] * n
        v[free] = 1
        for i, c in enumerate(pivots):
            v[c] = -rows[i][free] % p
        basis.append(v)
    return basis

def rank_into_quotient(vectors, sub):
    """Rank of the span of vectors in the quotient by span(sub)."""
    return rank(vectors + sub) - rank(sub)

def random_matrix(m, n, rng, p=P, sparse=0.5):
    """A random m x n matrix over F_p with many zeros (to get kernels)."""
    return [[rng.randrange(p) if rng.random() > sparse else 0
             for _ in range(n)] for _ in range(m)]

def snake_data(a, c, a2, c2, rng):
    """Random diagram; returns dims and ranks of the six-term sequence."""
    alpha = random_matrix(a2, a, rng)
    gamma = random_matrix(c2, c, rng)
    x = random_matrix(a2, c, rng, sparse=0.2)
    beta = ([alpha[i] + x[i] for i in range(a2)]
            + [[0] * a + gamma[j] for j in range(c2)])
    ka, kb, kc = (kernel_basis(alpha, a), kernel_basis(beta, a + c),
                  kernel_basis(gamma, c))
    im_a, im_b, im_c = columns(alpha), columns(beta), columns(gamma)
    f_ka = [v + [0] * c for v in ka]               # ker a -> ker b
    g_kb = [v[a:] for v in kb]                     # ker b -> ker c
    delta = [apply(x, w) for w in kc]              # ker c -> coker a
    fp = [e + [0] * c2 for e in identity(a2)]      # f' : A' -> B'
    gp = [e[a2:] for e in identity(a2 + c2)]       # g' : B' -> C'
    dims = [len(ka), len(kb), len(kc), a2 - rank(im_a),
            a2 + c2 - rank(im_b), c2 - rank(im_c)]
    ranks = [rank(f_ka), rank(g_kb), rank_into_quotient(delta, im_a),
             rank_into_quotient(fp, im_b), rank_into_quotient(gp, im_c)]
    zero_composites = (
        rank_into_quotient([apply(x, v[a:]) for v in kb], im_a) == 0
        and rank_into_quotient([apply(x, w) + [0] * c2 for w in kc], im_b)
        == 0)
    return dims, ranks, zero_composites

def is_exact(dims, ranks):
    """Exactness of 0 -> V1 -> ... -> V6 -> 0 from dims and map ranks."""
    incoming = [0] + ranks          # rank of the map into each term
    outgoing = ranks + [0]          # rank of the map out of each term
    return all(i == d - o for i, d, o in zip(incoming, dims, outgoing))

def main():
    rng = random.Random(2025)
    trials, exact, nonzero_delta, euler_ok = 300, 0, 0, 0
    example = None
    for t in range(trials):
        a, c, a2, c2 = (rng.randint(1, 4) for _ in range(4))
        dims, ranks, zero = snake_data(a, c, a2, c2, rng)
        exact += is_exact(dims, ranks) and zero
        nonzero_delta += ranks[2] > 0
        alt = dims[0] - dims[1] + dims[2] - dims[3] + dims[4] - dims[5]
        euler_ok += alt == 0
        if example is None and ranks[2] > 0 and min(dims) > 0:
            example = (a, c, a2, c2, dims, ranks)
    print(f"Snake lemma over F_{P}: {trials} random morphisms of short exact "
          "sequences")
    print(f"  six-term sequence exact at all six terms: {exact} / {trials}")
    print(f"  connecting map delta nonzero:             {nonzero_delta}")
    print(f"  alternating sum of dimensions equal to 0: {euler_ok} / {trials}")
    if example:
        a, c, a2, c2, dims, ranks = example
        names = ["ker a", "ker b", "ker c", "coker a", "coker b", "coker c"]
        print(f"Example with dim A, C, A', C' = {a}, {c}, {a2}, {c2}:")
        print("  dimensions: " + ", ".join(
            f"{n} = {d}" for n, d in zip(names[:3], dims[:3])))
        print("              " + ", ".join(
            f"{n} = {d}" for n, d in zip(names[3:], dims[3:])))
        print(f"  ranks of the five maps (delta in the middle): {ranks}")
    print()
    print("Epistemic status: exact arithmetic in F_5, so each check is a true "
          "statement\nabout the diagram tested; 300 random diagrams are "
          "evidence, not proof. The\nsnake lemma itself is proved in the "
          "lesson (and in Mathlib) for every\nabelian category.")

if __name__ == "__main__":
    main()

"""Where free abelian groups fail to be abelian: torsion in cokernels.

In the category F of finitely generated free abelian groups, a map
phi : Z^m -> Z^n has a kernel (subgroups of free groups are free), but its
cokernel in F is Z^n modulo the *saturation* of im(phi): the torsion of the
abelian-group cokernel Z^n / im(phi) is invisible. The canonical map
coimage -> image is then the inclusion im(phi) -> sat(im(phi)), of index
d_1 d_2 ... d_r (the invariant factors of phi), and it is an isomorphism
exactly when every invariant factor is 1. We compute invariant factors from
determinantal divisors (gcds of minors) for every 2 x 2 integer matrix with
entries in -2..2, and count how often F's coimage -> image map fails to be an
isomorphism - i.e. how often the abelian-category axiom fails.
"""

import itertools
import math

def det(mat):
    """Determinant of a small square integer matrix (Leibniz formula)."""
    n = len(mat)
    total = 0
    for perm in itertools.permutations(range(n)):
        sign = 1
        for i in range(n):
            for j in range(i + 1, n):
                if perm[i] > perm[j]:
                    sign = -sign
        term = sign
        for i in range(n):
            term *= mat[i][perm[i]]
        total += term
    return total

def determinantal_divisors(mat):
    """[D_1, D_2, ...]: D_k = gcd of all k x k minors (0 if all vanish)."""
    m, n = len(mat), len(mat[0])
    out = []
    for k in range(1, min(m, n) + 1):
        g = 0
        for rows in itertools.combinations(range(m), k):
            for cols in itertools.combinations(range(n), k):
                g = math.gcd(g, det([[mat[r][c] for c in cols] for r in rows]))
        out.append(g)
    return out

def invariant_factors(mat):
    """Nonzero invariant factors d_1 | d_2 | ... of an integer matrix."""
    divisors = [1] + [d for d in determinantal_divisors(mat) if d != 0]
    return [divisors[k] // divisors[k - 1] for k in range(1, len(divisors))]

def cokernels(mat):
    """(Ab-cokernel description, F-cokernel rank) for mat : Z^cols -> Z^rows.

    Ab: list of torsion orders d > 1 plus a free rank; F: the free rank only.
    """
    factors = invariant_factors(mat)
    free_rank = len(mat) - len(factors)
    torsion = [d for d in factors if d > 1]
    return (torsion, free_rank), free_rank

def coimage_to_image_index(mat):
    """Index of im(phi) in its saturation: the product of invariant factors."""
    return math.prod(invariant_factors(mat))

def describe(torsion, free_rank):
    """A readable name such as 'Z/2 + Z' for an abelian group."""
    parts = [f"Z/{d}" for d in torsion] + ["Z"] * free_rank
    return " + ".join(parts) if parts else "0"

def main():
    print("Multiplication by k on Z, as a map of free abelian groups:")
    for k in [1, 2, 3, 0]:
        (tors, free), f_rank = cokernels([[k]])
        print(f"  k = {k}: cokernel in Ab = {describe(tors, free):<4}"
              f" cokernel in F = {describe([], f_rank):<2}"
              f" coim -> im index = {coimage_to_image_index([[k]])}")
    entries = range(-2, 3)
    mats = [[[a, b], [c, d]] for a, b, c, d in itertools.product(entries,
                                                                 repeat=4)]
    by_index = {}
    for mat in mats:
        idx = coimage_to_image_index(mat)
        by_index[idx] = by_index.get(idx, 0) + 1
    bad = sum(n for idx, n in by_index.items() if idx > 1)
    print(f"All {len(mats)} integer 2 x 2 matrices with entries in -2..2:")
    for idx in sorted(by_index):
        print(f"  index [sat(im) : im] = {idx}: {by_index[idx]:3d} matrices")
    print(f"  coimage -> image not an isomorphism in F: {bad} matrices")
    example = [[2, 0], [0, 2]]
    (tors, free), _ = cokernels(example)
    print(f"  e.g. 2 * identity: cokernel in Ab = {describe(tors, free)}, "
          "cokernel in F = 0:")
    print("  monic and epic in F, but not invertible")
    print()
    print("Epistemic status: exact integer arithmetic, exhaustive over the "
          "625\nmatrices shown. It illustrates the theorem that F is additive "
          "but not\nabelian; the proof (and the Lean file) handle the general "
          "case.")

if __name__ == "__main__":
    main()