Smith Normal Form (pedagogical purposes, not efficiency)

SmithNormalForm.py

import numpy as np
from numpy.random import default_rng
rng=default_rng()

def rA(n):
    """make a random nxn integer matrix with entries in the range -30,30"""
    A=rng.integers(-30,30,size=n*n).reshape(n,n)
    return A

def swapcol(A,i,j):
    """interchange columns i and j in A"""
    B=A.copy()
    B[:,[i,j]]=B[:,[j,i]]
    return B

def swaprow(A,i,j):
    """interchange rows i and j in A"""
    B=A.copy()
    B[[i,j],:]=B[[j,i],:]
    return B

def elemcol(A,i,j,x):
    """elementary column operation: replace column i by column i - x* column j"""
    B=A.copy()
    B[:,i]=B[:,i]-x*B[:,j]
    return B
    
def elemrow(A,i,j,x):
    """elementary row operation: replace row i by row i -x * row j"""
    B=A.copy()
    B[i,:]=B[i,:]-x*B[j,:]
    return B

def reducecol(A,r,i,j):
    """do an elementary column operation reducing the entry in row r of column i mod row r of column j"""
    q=A[r,i]//A[r,j] 
    return elemcol(A,i,j,q)

def reducerow(A,r,i,j):
    """do an elementary row operation reducing the entry in column r of row i mod column r of row j"""
    q=A[i,r]//A[j,r]
    return elemrow(A,i,j,q)

def rowpivot(A,i,r):
    """swap columns so that the smallest nonzero entry in row i is in column r"""
    B=A.copy()
    R=np.abs(A[i,:])
    j=np.where(R==0,np.inf,R).argmin()
    if A[i,j]<0:
        B[i,:]=-B[i,:]
    return swapcol(B,j,r)

def colpivot(A,i,r):
    """swap rows so that the smallest nonzero intry in column i is in row r"""
    B=A.copy()
    R=np.abs(A[:,i])
    j=np.where(R==0,np.inf,R).argmin()
    if A[j,i]<0:
        B[:,i]=-B[:,i]
    return swaprow(B,j,r)

def smith(A):
    """put A in smith normal form"""
    B=A.copy()
    n=A.shape[0]
    for r in range(n-1):
        while True:
            bflag=True
            while np.max(np.abs(B[r,(r+1):]))>0:
                bflag=False
                B=rowpivot(B,r,r)
                for i in range(r+1,n):
                    B=reducecol(B,r,i,r)
            while np.max(np.abs(B[(r+1):,r]))>0:
                bflag=False
                B=colpivot(B,r,r)
                for i in range(r+1,n):
                    B=reducerow(B,r,i,r)
            if bflag:
                break
    B[n-1,n-1]=np.abs(B[n-1,n-1])
    return B