Smith normal form.

smith.py

import numpy as np

def smith_normal_form(matrix):
    A = np.array(matrix)
    m, n = A.shape
    P = np.eye(m)
    Q = np.eye(n)

    i = 0
    j = 0

    while i < m and j < n:
        # Find the pivot element
        pivot = A[i:, j:].argmin()
        pivot_i, pivot_j = divmod(pivot, n - j)
        pivot_i += i
        pivot_j += j
        pivot_value = A[pivot_i, pivot_j]

        # Swap rows
        if pivot_i != i:
            P[[i, pivot_i]] = P[[pivot_i, i]]
            A[[i, pivot_i]] = A[[pivot_i, i]]

        # Swap columns
        if pivot_j != j:
            Q[:, [j, pivot_j]] = Q[:, [pivot_j, j]]
            A[:, [j, pivot_j]] = A[:, [pivot_j, j]]

        # Perform row operations to make the pivot element the only non-zero element in its row and column
        for row in range(i + 1, m):
            factor = A[row, j] // pivot_value
            A[row, j:] -= factor * A[i, j:]
            P[row] -= factor * P[i]

        for col in range(j + 1, n):
            factor = A[i, col] // pivot_value
            A[i:, col] -= factor * A[i, j:]
            Q[:, col] -= factor * Q[:, j]

        i += 1
        j += 1

    return P, A, Q