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);print
show(transf(A,o));print
print(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