# Has Mathematica a function to compute the Smith Normal Form?

The Smith normal form is a matrix that can be calculated for any matrix (not necessarily square) with integer entries. See Wikipedia for a more elaborate description. Has Mathematica a function to compute the Smith Normal Form ? If not, is there an easy way to compute it in Mathematica?

-

I don't think Mathematica has a built-in function to compute the Smith Normal Form, but here are a couple packages by David Jabon that can do this for you: Smith Normal Form

Usage

Needs@"IntegerSmithNormalForm"


Here is an integer matrix mat

mat = {{1, 2, 3}, {-2, 3, 1}, {3, 2, 1}}


We compute it's SmithForm

smith = SmithForm[mat]


One can even compute the ExtendedSmithForm

(* "ExtendedSmithForm[A] gives, for an integral
matrix A, {D, {P,Q}} where
D is in Smith normal form and P and Q are
matrices such that P A Q = D. "  *)

extsmith = ExtendedSmithForm[mat]


Gives:

{{{1, 0, 0}, {0, 1, 0}, {0, 0, 28}},
{{{1, 0, 0}, {4, -1, -2}, {13, -4, -7}}, {{1, -2, 15}, {0, 1, -9}, {0, 0, 1}}}}


We can check that indeed P A Q == D

extsmith[[2, 1]].mat.extsmith[[2, 2]] == smith


True

-
Coded in 1994 and it still works in Mathematica 9, that's how it should be. –  ndroock1 Jan 22 at 18:59
@ndroock1. I agree. –  RunnyKine Jan 22 at 19:01

Other than indenting the code blocks, what follows comes verbatim from an old file of mine. I have not tested it in the current Mathematica, though I have no reason to think it might fail.

diagonalQ[mat_?MatrixQ] := With[
{posns=Flatten[Map[Position[#, _?(#!=0&)]&, mat]]},
Length[Union[posns]] == Length[posns]
]

diagonalize[mat_?MatrixQ] := Module[
{hnf=mat, umat=IdentityMatrix[Length[mat]],
vmat=IdentityMatrix[Length[mat[[1]]]], tmpu, tmpv},
While[Not[diagonalQ[hnf]],
{tmpu, hnf} = HermiteDecomposition[hnf];
umat = tmpu . umat;
{tmpv, hnf} = HermiteDecomposition[Transpose[hnf]];
vmat = vmat . Transpose[tmpv];
hnf = Transpose[hnf];
];
{umat, hnf, vmat}
]

divides[a_, b_] := GCD[a,b]===a

smithNormalForm[mat_?MatrixQ] := Module[
{uu, dd, vv, diags, gcd, col=0, dim, tmpu, tmpv},
{uu, dd, vv} = diagonalize[mat];
diags = Select[Flatten[dd], #!=0&];
dim = Length[diags];
While[col+1<dim,
col++;
If [divides[diags[[col]],GCD[Apply[Sequence,Drop[diags,col]]]],
Continue[]];
vv = Transpose[vv];
Do [
dd[[j,col]] = diags[[j]];
vv[[col]] += vv[[j]],
{j,col+1,dim}];
vv = Transpose[vv];
{tmpu, dd, tmpv} = diagonalize[dd];
uu = tmpu  . uu;
vv = vv . tmpv;
diags = Select[Flatten[dd], #!=0&];
];
{uu, dd, vv}
]


For examples:

SeedRandom[1111];
m = 4; n = 6;
uppermat[m_] :=
Table[If[i==j,1,If[i>j,0,Random[Integer,{-10,10}]]], {i,m}, {j,m}]
umat[m_] := uppermat[m] . Transpose[uppermat[m]]
uu = umat[m];
vv = umat[n];
hh = {{1,0,0,0,0,0}, {0,3,0,0,0,0}, {0,0,0,18,0,0}, {0,0,0,0,0,90}};
testmat = uu . hh . vv;


To change this just rig hh so it is appropriate for the dimensions.
To get nontrivial examples you might want to have some rows and/or
columns zeroed, and you'd want nontrivial gcds for the nonzero
elements.

{newu, newh, newv} = smithNormalForm[testmat]
Max[Abs[newu . testmat . newv - newh]] (* want zero *)

{newu, newh, newv} = smithNormalForm[Transpose[testmat]]
Max[Abs[newu . Transpose[testmat] . newv - newh]] (* want zero *)
`
-
I really wonder why Mma does not have this builtin. –  Mariano Suárez-Alvarez Jun 28 at 19:18
@Mariano Suárez-Alvarez About all I can say is I wonder about it myself. –  Daniel Lichtblau Jun 28 at 20:05