Has Mathematica a function to compute the Smith Normal Form?

Question

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?

Answer 1

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

Answer 2

Just to separate this from a package-based answer. In Mathematica 10.2, you can now do this with the built-in function SmithDecomposition. So using the same matrix from my previous answer:

mat = {{1, 2, 3}, {-2, 3, 1}, {3, 2, 1}};
MatrixForm /@ SmithDecomposition[mat]

Where the second element is the Smith normal form.

Answer 3

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 *)