# 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. Mar 12, 2014 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). May 14, 2020 at 15:14
• I guess Daniel is suggesting something like SmithDecomposition[ArrayFlatten[{{mat, 2 IdentityMatrix[nrows]}, {2 IdentityMatrix[ncols], 0}}]]. May 14, 2020 at 15:30
• What @J.M... said. If I knew how to use ArrayFlatten. May 14, 2020 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]];
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).

• Daniel, your knowledge of "what has already been done" seems to be amazing. Mar 12, 2014 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.) Mar 12, 2014 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? Mar 12, 2014 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. Mar 12, 2014 at 16:37