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.

I would like to compute the Smith Normal form of a matrix with coefficients in GF(p). In particular, I am interested in GF(2). I have used the Smith normal form packages for integer and polynomial matrices before. Can these be adapted to work for matrices with finite field coefficients? Is there another way to compute the Smith normal forms mod p in Mathematica?

share|improve this question
    
A small example input would be helpful. –  Daniel Lichtblau Mar 12 at 15:24

1 Answer 1

Here is some code based on functionality in the development version of Mathematica.

moduleGroebnerBasis[polys_, p_, vars_, cvars_, opts___] := 
 Catch[Module[{newpols, rels, len = Length[cvars], gb, j, k, rul}, 
   If[! FreeQ[polys, $Failed], Throw[$Failed]];
   rels = Flatten[Table[cvars[[j]]*cvars[[k]], {j, len}, {k, j, len}]];
   newpols = Join[polys, rels];
   Quiet[gb = 
     GroebnerBasis[newpols, Join[cvars, vars], opts, Modulus -> p]];
   If[Head[gb] === GroebnerBasis, Throw[$Failed]];
   rul = Map[(# -> {}) &, rels];
   gb = Flatten[gb /. rul];
   Collect[gb, cvars]]]

groebnerHNF[omat_?MatrixQ, p_, var_] := 
 Catch[Module[{mat, nr, nc, v, newvars, generators, mgb, 
    res}, {nr, nc} = Dimensions[omat];
   mat = Join[omat, IdentityMatrix[nr], 2];
   newvars = Array[v, nr + nc];
   generators = mat.newvars;
   mgb = moduleGroebnerBasis[generators, p, {var}, newvars];
   If[! FreeQ[mgb, $Failed] || Length[mgb] =!= nr, Throw[$Failed]];
   res = Outer[D, Reverse[mgb], newvars];
   {res[[All, 1 ;; nc]], res[[All, nc + 1 ;; -1]]}]]

diagonalMatrixQ[mat_?MatrixQ] := 
 Catch[Do[If[i == j, Continue[], 
     If[mat[[i, j]] =!= 0, Throw[False]]];, {i, Length[mat]}, {j, 
    Length[mat[[1]]]}];
  Throw[True]]

diagonalize[mat_, p_, var_] := 
 Catch[Module[{hnf = mat, nr = Length[mat], nc = Length[mat[[1]]], 
    umat, vmat, tmpu, tmpv, approx = Precision[mat] =!= Infinity}, 
   umat = IdentityMatrix[nr];
   vmat = IdentityMatrix[nc];
   While[Not[diagonalMatrixQ[hnf]], hnf = groebnerHNF[hnf, p, var];
    If[! FreeQ[hnf, $Failed], Throw[$Failed]];
    {hnf, tmpu} = hnf;
    umat = PolynomialMod[Dot[tmpu, umat], p];
    hnf = groebnerHNF[Transpose[hnf], p, var];
    If[! FreeQ[hnf, $Failed], Throw[$Failed]];
    {hnf, tmpv} = hnf;
    vmat = PolynomialMod[Dot[vmat, Transpose[tmpv]], p];
    hnf = Transpose[hnf];];
   {umat, hnf, vmat}]]

dividesQ[p1_, p2_, p_, var_] := Catch[Module[{quo, rem},
   If[FreeQ[{p1, p2}, var], Throw[True]];
   {quo, rem} = PolynomialQuotientRemainder[p2, p1, var, Modulus -> p];
   rem === 0]]

smithDecompositionGroebner[mat_, p_, var_] := 
 Catch[Module[{snf, uu, dd, vv, diags, gcd, col = 0, dim, tmpu, tmpv},
    snf = diagonalize[mat, p, var];
   If[! FreeQ[snf, $Failed], Throw[$Failed]];
   {uu, dd, vv} = snf;
   diags = Select[Flatten[dd], Not[TrueQ[# == 0]] &];
   dim = Length[diags];
   While[col + 1 < dim, col++;
    If[dividesQ[diags[[col]], 
      PolynomialGCD[Apply[Sequence, Drop[diags, col]], Modulus -> p], 
      p, var], Continue[]];
    vv = Transpose[vv];
    Do[dd[[j, col]] = diags[[j]];
     vv[[col]] = PolynomialMod[vv[[col]] + vv[[j]], p], {j, col + 1, 
      dim}];
    vv = Transpose[vv];
    snf = diagonalize[dd, p, var];
    If[! FreeQ[snf, $Failed], Throw[$Failed]];
    {tmpu, dd, tmpv} = snf;
    uu = PolynomialMod[tmpu.uu, p];
    vv = PolynomialMod[vv.tmpv, p];
    diags = Select[Flatten[dd], Not[TrueQ[# == 0]] &];];
   {uu, dd, vv}]]

Brief example:

t = {{(1 + s) (3 + s), 0, 0}, {(1 + s) (3 + s), 1 + s, 0}, {0,
     1 + s, (1 + s) (4 + s)}};

smithDecompositionGroebner[t, 13, s]

(* Out[80]= {{{0, 12, 1}, {12, 2, 12}, {12, 10 + 12 s, 3 + s}}, {{1 + s, 
   0, 0}, {0, 1 + s, 0}, {0, 0, 12 + 6 s + 8 s^2 + s^3}}, {{1, 0, 
   9 + 12 s}, {1, 1, 0}, {1, 0, 10 + 12 s}}} *)

I perfectly reasonable question might be "Does it work correctly?" My response would be "Probably". The basic approach to the Hermite decomposition is described here (I seem to be referring to this quite a bit today).

share|improve this answer
    
Daniel, your knowledge of "what has already been done" seems to be amazing. –  Jacob Akkerboom Mar 12 at 16:03
    
@Jacob Akkerboom It was easy enough in this case, since it involved (mostly) code I had written. (I'm not a Control Theorist, but sometimes I portray one in my spare time.) –  Daniel Lichtblau Mar 12 at 16:17
    
I was just wondering why there is still no build-in function in reference.wolfram.com/language/guide/MatrixPredicates.html for the function you wrote above diagonalMatrixQ? Should not this predict be also a build-in function in Mathematica? –  Nasser Mar 12 at 16:32
    
@Nasser Hard for me to say whether it warrants built-in status. I can say it's not a predicate I've used often myself. But maybe others do, and frequently write their own. –  Daniel Lichtblau Mar 12 at 16:37

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.