# Computing Smith Normal form of a matrix with mod p coefficients

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?

-
A small example input would be helpful. – Daniel Lichtblau Mar 12 '14 at 15:24

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

-
Daniel, your knowledge of "what has already been done" seems to be amazing. – Jacob Akkerboom Mar 12 '14 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 '14 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 '14 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 '14 at 16:37