def showop(o): l="" for op in o: if op[0]==1:l=l+"C"+str(op[1])+" + "+str(op[3])+"C"+str(op[2])+" → C"+str(op[1])+"\n" elif op[0]==2:l=l+'L'+str(op[1])+' + '+str(op[3])+'*L'+str(op[2])+' → L'+str(op[1])+"\n" elif op[0]==3:l=l+'C'+str(op[1])+" ↔ C"+str(op[2])+"\n" elif op[0]==4:l=l+'L'+str(op[1])+" ↔ L"+str(op[2])+"\n" return(l)
import copy def transf(A,o): B=copy.deepcopy(A) for op in o: if op[0]==1:B[:,op[1]]=B[:,op[1]]+op[3]*B[:,op[2]] elif op[0]==2:B[op[1],:]=B[op[1],:]+op[3]*B[op[2],:] elif op[0]==3: C=B[:,op[1]];B[:,op[1]]=B[:,op[2]];B[:,op[2]]=C elif op[0]==4: L=B[op[1],:];B[op[1],:]=B[op[2],:];B[op[2],:]=L return(B)
def fsmith(A): p=A.nrows();q=A.ncols() L=[abs(A[i][j]) for i in range(p) for j in range(q)] if p*q==0 or max(L)==0: return([]) else: a=min(L[i] for i in range(p*q) if L[i]!=0) k=L.index(a);i0=k//q;j0=k%q;a=A[i0][j0] l=[i for i in range(p*q) if L[i]%a!=0] if l==[]: o=[[1,j,i0,-A[i0][j]/a] for j in [0..q-1] if j!=j0]+[[2,i,j0,-A[i][j0]/a] for i in [0..p-1] if i!=i0] if j0!=0:o=o+[[3,0,j0]] if i0!=0:o=o+[[4,0,i0]] return(o+[[op[0],1+op[1],1+op[2]]+op[3:] for op in fsmith(transf(A,o).submatrix(1,1))]) else: k=min(l);i1=k//q;j1=k%q;b=A[i1][j1] if A[i0][j1]%a != 0: o=[[1,j1,j0,-(A[i0][j1]//a)]] elif A[i1][j0]%a != 0: o=[[2,i1,i0,-(A[i1][j0]//a)]] else: o=[[2,i1,i0,-(A[i1][j0]//a)],[2,i0,i1,1],[1,j1,j0,-(((1-A[i1][j0]//a)*A[i0][j1]+b)//a)]] return(o+fsmith(transf(A,o)))
A=matrix(2,3,[4,2,4,3,6,6])o=fsmith(A)show(A);printshow(transf(A,o));printprint(showop(o))
$\displaystyle \left(\begin{array}{rrr}
4 & 2 & 4 \\
3 & 6 & 6
\end{array}\right)$
$\displaystyle \left(\begin{array}{rrr}
1 & 0 & 0 \\
0 & 6 & 0
\end{array}\right)$
L1 + -3*L0 → L1
L0 + 1*L1 → L0
C0 + 3C1 → C0
C1 + -2C0 → C1
C2 + 2C0 → C2
L1 + 9*L0 → L1
C2 + 2C1 → C2
C1 + -1C2 → C1
C2 + -2C1 → C2