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

2 Answers 2

up vote 7 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 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]

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

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

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.