# Computing Smith normal form of a matrix with $\bmod 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 $\bmod p$ in Mathematica?

• A small example input would be helpful. – Daniel Lichtblau Mar 12 '14 at 15:24
• Since this has bubbled back to the surface, I think today what I might try is this. Augment beneath the matrix with 2*IdentityMatrix[ncols], to the right with 2*IdentityMatrix[nrows], and the bottom right with zeros of dimension ncols x nrows. Then do SmithDecomposition. Last, remove excess (this might take some figuring out though). – Daniel Lichtblau May 14 at 15:14
• I guess Daniel is suggesting something like SmithDecomposition[ArrayFlatten[{{mat, 2 IdentityMatrix[nrows]}, {2 IdentityMatrix[ncols], 0}}]]. – J. M.'s technical difficulties May 14 at 15:30
• What @J.M... said. If I knew how to use ArrayFlatten. – Daniel Lichtblau May 14 at 17:58

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]];
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[]]}];
Throw[True]]

diagonalize[mat_, p_, var_] :=
Catch[Module[{hnf = mat, nr = Length[mat], nc = Length[mat[]],
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= {{{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).

• 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