# Find the order $m$ of a matrix ${\bf A}$ such that ${\bf A}^m= {\bf 1}$

I have a square matrix ${\bf A}$ defined over the field $\mathbb Z_2$ and I want to find its order such that ${\bf A}^m=1$. I tried using

Solve[MatrixPower[A,m] == IdentityMatrix[Length[A]], m, Modulus-> 2]


but this seems to be taking forever. How should I find the order of the matrix as defined? For example, A = {{1, 0}, {1, 1}} has order $2$.

• Could you show a small typical matrix so that people might experiment, confirm what you are seeing and try possible solutions before finding that they are doing something very different from you? – Bill May 8 '18 at 20:46
• Maybe something like: First@Position[ Table[Mod[MatrixPower[A, n], 2] == IdentityMatrix[2], {n, 4}], True]? Can be improved by stopping computations when True is detected. – anderstood May 8 '18 at 21:06
• Or: B = A; n = 0; While[(B != IdentityMatrix[2]) && (n < 100), B = Modulus[A.B, 2]; n++]; n – anderstood May 8 '18 at 21:10
• m = 1; While[m<=100&& Mod[MatrixPower[A, m++], 2] =!= IdentityMatrix[Length[A]]]; m - 1? – kglr May 8 '18 at 21:13
• also mmax = 1000; Block[{m = 1, B = A}, While[m++ <= mmax && (B = Mod[B.A, 2]) =!= IdentityMatrix[Length[A]]]; m] – kglr May 9 '18 at 1:01

Using the example in the OP,

A = {{1, 0}, {1, 1}};


we compare the different options mentioned in the comments (plus one that is mine):

Table[If[Mod[MatrixPower[A, n], 2] == IdentityMatrix[Length@A], Return[n, Table]], {n, 1000}] // RepeatedTiming
(* {0.00007, 2}, by anderstood *)

Block[{m = 1}, While[m <= 100 && Mod[MatrixPower[A, m++], 2] =!= IdentityMatrix[Length[A]]]; m - 1] // RepeatedTiming
(* {0.00003061, 2}, by kglr *)

Length[NestWhileList[Mod[#.A, 2] &, A, Mod[#, 2] =!= IdentityMatrix[Length[A]] &]] // RepeatedTiming
(* {9.8*10^-6, 2}, by me *)

Block[{mmax = 1000, m = 1, B = A}, While[m++ <= mmax && (B = Mod[B.A, 2]) =!= IdentityMatrix[Length[A]]]; m] // RepeatedTiming
(* {5.5*10^-6, 2}, by kglr *)


For less trivial examples, consider

SeedRandom[2];
A = RandomInteger[{0, 1}, {5, 5}];
(* {0.001, 12}, by anderstood *)
(* {0.00088, 12}, by kglr *)
(* {0.000064, 12}, by me *)
(* {0.000046, 12}, by kglr *)

SeedRandom[2];
A = RandomInteger[{0, 1}, {10, 10}];
(* {0.04, 60}, by anderstood *)
(* {0.037, 60}, by kglr *)
(* {0.00051, 60}, by me *)
(* {0.00037, 60}, by kglr *)

• Changing Modulus to Mod gives basically the same answer as kglr number 2, and the timing is very similar: Block[{}, B = A; n = 0; While[(B =!= IdentityMatrix[5]) && (n < 1000), B = Mod[A.B, 2]; n++]; n + 1]. And here is an updated version of my first comment with a break: Table[If[Mod[MatrixPower[A, n], 2] == IdentityMatrix[Length@A], Return[n, Table]], {n, 1000}] // RepeatedTiming. – anderstood May 9 '18 at 15:20

## The Code

I solve the more generic question of finding the power of a matrix under any mod-z with the following command:

Options[solver] = {matrixForm -> False};
solver[input_, mod_, iteration_, opts : OptionsPattern[]] :=
Module[{a, exp, exp2},
a /: Power[a, n_] := a;
a /: Times[a, n_] /; (n > mod - 1) := Mod[n, mod] a;
exp = Simplify[a input];
exp2 = NestWhileList[{#[[1]] + 1, exp.#[[2]]} &, {1, exp},
And[(! Simplify[#[[2]] === a IdentityMatrix[Length@input]]), #[[
1]] < iteration] &];
If[OptionValue[matrixForm] == True,
Grid[{#[[1]], #[[2]] // MatrixForm} & /@ (exp2 /. a -> 1)],
exp2 /. a -> 1
]
];


The code basically multiplies the given matrix with a scooping variable which is constrained such that its coefficients are in modulo-chosen number (second argument). Then we calculate its matrix powers till we find the identity matrix or we reach the chosen iteration number (third argument). We also have the option to make the output more visual.

## Check of OP's matrix

solver[{{1, 0}, {1, 1}}, 2, 100] // Timing


{0.000244, {{1, {{1, 0}, {1, 1}}}, {2, {{1, 0}, {0, 1}}}}}

where we get a list of matrices $\{A,A^2,A^3,\cdots\}$. We can make it more visual as follows:

solver[{{1, 0}, {1, 1}}, 2, 100, matrixForm -> True]


$$\begin{array}{cc} 1 & \left( \begin{array}{cc} 1 & 0 \\ 1 & 1 \\ \end{array} \right) \\ 2 & \left( \begin{array}{cc} 1 & 0 \\ 0 & 1 \\ \end{array} \right) \\ \end{array}$$

We can also check modulo other numbers:

solver[{{1, 0}, {1, 1}}, 3, 100, matrixForm -> True] // Timing


$$\left\{0.000352, \begin{array}{cc} 1 & \left( \begin{array}{cc} 1 & 0 \\ 1 & 1 \\ \end{array} \right) \\ 2 & \left( \begin{array}{cc} 1 & 0 \\ 2 & 1 \\ \end{array} \right) \\ 3 & \left( \begin{array}{cc} 1 & 0 \\ 0 & 1 \\ \end{array} \right) \\ \end{array} \right\}$$

## Another example

A = RandomInteger[{0, 1}, {5, 5}]


{{1, 0, 0, 0, 1}, {0, 1, 0, 1, 0}, {0, 1, 0, 0, 1}, {0, 1, 1, 1, 0}, {0, 0, 0, 1, 0}}

solver[A, 2, 100, matrixForm -> True] // Timing


$$\left\{0.001191, \begin{array}{cc} 1 & \left( \begin{array}{ccccc} 1 & 0 & 0 & 0 & 1 \\ 0 & 1 & 0 & 1 & 0 \\ 0 & 1 & 0 & 0 & 1 \\ 0 & 1 & 1 & 1 & 0 \\ 0 & 0 & 0 & 1 & 0 \\ \end{array} \right) \\ 2 & \left( \begin{array}{ccccc} 1 & 0 & 0 & 1 & 1 \\ 0 & 0 & 1 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 \\ 0 & 1 & 1 & 0 & 1 \\ 0 & 1 & 1 & 1 & 0 \\ \end{array} \right) \\ 3 & \left( \begin{array}{ccccc} 1 & 1 & 1 & 0 & 1 \\ 0 & 1 & 0 & 0 & 1 \\ 0 & 1 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 1 \\ 0 & 1 & 1 & 0 & 1 \\ \end{array} \right) \\ 4 & \left( \begin{array}{ccccc} 1 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 1 \\ \end{array} \right) \\ \end{array} \right\}$$

For the same matrix, we can find the necessary power in module-3 as well:

solver[A, 3, 100] // Timing


{0.006589, {{1, {{1, 0, 0, 0, 1}, {0, 1, 0, 1, 0}, {0, 1, 0, 0, 1}, {0, 1, 1, 1, 0}, {0, 0, 0, 1, 0}}}, {2, {{1, 0, 0, 1, 1}, {0, 2, 1, 2, 0}, {0, 1, 0, 2, 0}, {0, 0, 1, 2, 1}, {0, 1, 1, 1, 0}}}, {3, {{1, 1, 1, 2, 1}, {0, 2, 2, 1, 1}, {0, 0, 2, 0, 0}, {0, 0, 2, 0, 1}, {0, 0, 1, 2, 1}}}, {4, {{1, 1, 2, 1, 2}, {0, 2, 1, 1, 2}, {0, 2, 0, 0, 2}, {0, 2, 0, 1, 2}, {0, 0, 2, 0, 1}}}, {5, {{1, 1, 1, 1, 0}, {0, 1, 1, 2, 1}, {0, 2, 0, 1, 0}, {0, 0, 1, 2, 0}, {0, 2, 0, 1, 2}}}, {6, {{1, 0, 1, 2, 2}, {0, 1, 2, 1, 1}, {0, 0, 1, 0, 0}, {0, 0, 2, 2, 1}, {0, 0, 1, 2, 0}}}, {7, {{1, 0, 2, 1, 2}, {0, 1, 1, 0, 2}, {0, 1, 0, 0, 1}, {0, 1, 2, 0, 2}, {0, 0, 2, 2, 1}}}, {8, {{1, 0, 1, 0, 0}, {0, 2, 0, 0, 1}, {0, 1, 0, 2, 0}, {0, 0, 0, 0, 2}, {0, 1, 2, 0, 2}}}, {9, {{1, 1, 0, 0, 2}, {0, 2, 0, 0, 0}, {0, 0, 2, 0, 0}, {0, 0, 0, 2, 0}, {0, 0, 0, 0, 2}}}, {10, {{1, 1, 0, 0, 1}, {0, 2, 0, 2, 0}, {0, 2, 0, 0, 2}, {0, 2, 2, 2, 0}, {0, 0, 0, 2, 0}}}, {11, {{1, 1, 0, 2, 1}, {0, 1, 2, 1, 0}, {0, 2, 0, 1, 0}, {0, 0, 2, 1, 2}, {0, 2, 2, 2, 0}}}, {12, {{1, 0, 2, 1, 1}, {0, 1, 1, 2, 2}, {0, 0, 1, 0, 0}, {0, 0, 1, 0, 2}, {0, 0, 2, 1, 2}}}, {13, {{1, 0, 1, 2, 0}, {0, 1, 2, 2, 1}, {0, 1, 0, 0, 1}, {0, 1, 0, 2, 1}, {0, 0, 1, 0, 2}}}, {14, {{1, 0, 2, 2, 2}, {0, 2, 2, 1, 2}, {0, 1, 0, 2, 0}, {0, 0, 2, 1, 0}, {0, 1, 0, 2, 1}}}, {15, {{1, 1, 2, 1, 0}, {0, 2, 1, 2, 2}, {0, 0, 2, 0, 0}, {0, 0, 1, 1, 2}, {0, 0, 2, 1, 0}}}, {16, {{1, 1, 1, 2, 0}, {0, 2, 2, 0, 1}, {0, 2, 0, 0, 2}, {0, 2, 1, 0, 1}, {0, 0, 1, 1, 2}}}, {17, {{1, 1, 2, 0, 2}, {0, 1, 0, 0, 2}, {0, 2, 0, 1, 0}, {0, 0, 0, 0, 1}, {0, 2, 1, 0, 1}}}, {18, {{1, 0, 0, 0, 0}, {0, 1, 0, 0, 0}, {0, 0, 1, 0, 0}, {0, 0, 0, 1, 0}, {0, 0, 0, 0, 1}}}}}

where we see that it is 18th power, and we get the result in around 6 ms, which is fast enough.

## Another example

A = RandomInteger[{0, 4}, {4, 4}]


{{2, 1, 1, 1}, {1, 2, 0, 2}, {4, 3, 2, 2}, {1, 4, 3, 2}}

With[{result = solver[A, 5, 100] // Timing}, {result[[1]],Last[result[[2]]]}]


$$\left\{0.02803,\left\{62,\left( \begin{array}{cccc} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ \end{array} \right)\right\}\right\}$$

We can compare the result with direct computation:

Mod[#, 5] & /@ Nest[#.A &, A, 61]


{{1, 0, 0, 0}, {0, 1, 0, 0}, {0, 0, 1, 0}, {0, 0, 0, 1}}

which proves that indeed 62nd power of the matrix A above is identity matrix for modulo-5 integers.

• For such short times, I think it's better to use RepeatedTiming instead of Timing (not that it's important to get an accurate timing here, just saying...) – AccidentalFourierTransform May 9 '18 at 18:30
• @AccidentalFourierTransform Thanks! – Soner May 10 '18 at 2:26

Another implementation of Soner's solver:

cf = Compile[{{A, _Integer, 2}, {B, _Integer, 2}, {mod, _Integer, 0}},
Mod[A.B, mod], Parallelization -> True]
solver[A_, mod_, maxIter_: 1000] :=
Block[{id = IdentityMatrix[Length@A], B = A, n = 0},
While[(B =!= id) && (n < 1000), B = cf[A, B, mod]; n++]; n + 1]


Examples:

A = {{2, 1, 1, 1}, {1, 2, 0, 2}, {4, 3, 2, 2}, {1, 4, 3, 2}};
solver[A, 5] // RepeatedTiming
(* {0.00022, 62} *)

A = {{1, 0, 0, 0, 1}, {0, 1, 0, 1, 0}, {0, 1, 0, 0, 1}, {0, 1, 1, 1, 0},
{0, 0, 0, 1, 0}};
solver[A, 2] // RepeatedTiming
(* {0.000018, 4} *)

solver[A, 3] // RepeatedTiming
(* {0.000074, 18} *)

• @AccidentalFourierTransform You're right. And with RepeatedTiming the time is greatly reduced; there must be some overhead. – anderstood May 9 '18 at 18:45

MatrixPowerMod from A Course in Computational Number Theory by Bressoud and Wagon, page 34, is quite fast.

MatrixPowerMod[a_, m_, p_] :=
Block[{b = a, d = IntegerDigits[m, 2]},
Do[
b = Mod[b.b, p];
If[d[[i]] == 1, b = Mod[b.a, p]],
{i, 2, Length[d]}];
b]


Then use some information from an answer to this question by loup blanc.

MatrixOrderMod[a_, p_] :=
Block[{d, s, k, ii = IdentityMatrix[Length[a]]},
d = Mod[Det[a], p];
If[d == 0, 0,
s = MultiplicativeOrder[d, p];
k = s;
While[MatrixPowerMod[a, k, p] =!= ii, k += s];
k]
]


MatrixOrderMod is a brute force search, as others have implemented, but is faster in general because the matrix exponent is incremented by s instead of by 1. The code also tests for zero determinant, which corresponds to no solution. This test eliminates the need for an explicit upper bound on the exponent.

• Nice. I thought about the determinant test too, but I wasn't sure: there is a solution if and only if $\det(A)\neq 0$? the only if direction is obvious, but the other one -- not quite so, at least to me. – AccidentalFourierTransform May 9 '18 at 23:31
• This can be compacted a bit: MatrixOrderMod[a_, p_] := Block[{k = 0, ii = IdentityMatrix[Length[a]], d, s}, d = Det[a, Modulus -> p]; If[d == 0, 0, s = MultiplicativeOrder[d, p]; While[AlgebraMatrixPowerMod[a, k += s, p] =!= ii]; k]]. – J. M.'s torpor Oct 1 '18 at 20:32