Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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?

share|improve this question
add comment

2 Answers

up vote 6 down vote accepted

I don't think Mathematica has a built-in function to compute the Smith Normal Form, but here are a couple packages 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]

Mathematica graphics

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

share|improve this answer
    
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
add comment

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 *)
share|improve this answer
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.