Computing Smith normal form of a matrix with $\bmod p$ coefficients

Question

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 $\bmod p$ in Mathematica?

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

A 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).