I will use $p$ for $13$ for short. Let $F$ be the field with $p=13$ elements, $F=\Bbb Z/13=\Bbb F_{13}$.
Let $R=\Bbb Z/13^2$ be the ring of integers modulo $p^2=13^2=169$. There is a canonical map $R\to F$, it takes $r\in R$ and goes modulo $p$, it is the map $\Bbb Z/p^2\Bbb Z\to\Bbb Z/p\Bbb Z$.
Now let us denote by $M$ the given matrix with coefficients in $R$, $C\in M_{2p\times 2p}(R)$.
We have then a projected version $S$ of it inside $M_{2p\times 2p}(F)$.
Let us first find the multiplicative order of $S$, since sage was educated only to do such computations
over a field. Denote by $n$ this order, so $S^n=I$ over $F$. This means that over the ring of matrices with
coefficients in $R$ we have an equality of the shape:
$$
M^n = I+pT\ .
$$
Then the binomial formula gives us immediately
$$
M^{pn}=(I+pT)^p=I+\binom p1 pT+p^2(\dots)=I\ .
$$
So the order $m$ of $M$ is either $n$ or $pn$.
Now we can go down by successively dividing $m$ by primes in $m$
as long as we get a diagonal matrix which is a multiple of the identity.
(If such a prime in $m$ leads to an alien form, we reject it in the following steps.)
There must be something strange with the definition of $M$ in the post.
The code for the story is as follows (and please format it for the human eye if possible next time, the eye is also thinking about the situation and a solution):
p = 13
F = Integers(p)
R = Integers(p^2)
A = matrix(
ZZ,
[
[ 20, 101, 52, 52, 166, 148, 46, 135, 96, 51, 73, 49, 128],
[166, 164, 159, 66, 123, 50, 144, 85, 29, 116, 22, 93, 10],
[158, 152, 90, 65, 20, 167, 27, 96, 109, 154, 127, 164, 76],
[120, 154, 132, 110, 22, 113, 115, 51, 25, 104, 108, 82, 33],
[ 43, 148, 131, 45, 81, 2, 164, 145, 117, 157, 4, 108, 61],
[134, 23, 151, 120, 151, 44, 30, 1, 76, 32, 60, 132, 165],
[121, 40, 83, 4, 56, 88, 3, 134, 100, 85, 88, 18, 3],
[ 23, 20, 20, 31, 66, 24, 41, 126, 47, 137, 33, 112, 49],
[143, 18, 44, 26, 89, 109, 118, 148, 35, 16, 35, 122, 150],
[144, 51, 47, 143, 109, 164, 52, 38, 92, 50, 98, 60, 104],
[ 70, 165, 89, 80, 28, 75, 19, 110, 101, 41, 155, 78, 67],
[123, 147, 54, 4, 60, 133, 49, 151, 30, 32, 157, 108, 82],
[ 85, 139, 50, 70, 124, 168, 87, 63, 13, 104, 58, 107, 113],
]
);
B = matrix.identity(ZZ, p)
C = 73*b
D = matrix.zero(ZZ, p, p)
Z = block_matrix([[A, C], [B, D]])
S = Z.base_extend(F)
M = Z.base_extend(R)
n = S.multiplicative_order()
print(f"S has order {n} = {n.factor()}")
And we obtain:
S has order 67575144 = 2^3 * 3 * 7 * 13 * 30941
It turns out that $M^n$ is not the identity matrix, it is instead:
$$
\scriptscriptstyle
\begin{bmatrix}
1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
91 & 1 & 156 & 143 & 78 & 130 & 156 & 13 & 39 & 156 & 0 & 117 & 0 & 13 & 39 & 117 & 156 & 78 & 143 & 65 & 117 & 156 & 52 & 78 & 65 & 0 \\
65 & 65 & 53 & 156 & 117 & 156 & 117 & 143 & 117 & 156 & 104 & 39 & 0 & 13 & 91 & 26 & 104 & 13 & 143 & 156 & 52 & 104 & 91 & 52 & 156 & 0 \\
117 & 156 & 0 & 1 & 130 & 52 & 52 & 143 & 130 & 104 & 26 & 143 & 0 & 143 & 26 & 91 & 0 & 65 & 104 & 117 & 156 & 39 & 65 & 143 & 26 & 0 \\
13 & 130 & 0 & 0 & 53 & 156 & 156 & 91 & 52 & 143 & 78 & 91 & 0 & 91 & 78 & 104 & 0 & 26 & 143 & 13 & 130 & 117 & 26 & 91 & 78 & 0 \\
91 & 52 & 65 & 0 & 26 & 14 & 143 & 130 & 104 & 117 & 156 & 13 & 0 & 117 & 26 & 143 & 104 & 91 & 91 & 0 & 104 & 52 & 117 & 143 & 65 & 0 \\
13 & 104 & 130 & 78 & 65 & 78 & 14 & 104 & 156 & 156 & 156 & 0 & 0 & 52 & 26 & 104 & 117 & 143 & 91 & 65 & 130 & 52 & 156 & 130 & 104 & 0 \\
104 & 52 & 39 & 130 & 156 & 0 & 156 & 131 & 117 & 78 & 130 & 130 & 0 & 0 & 156 & 65 & 91 & 156 & 52 & 0 & 156 & 130 & 39 & 104 & 78 & 0 \\
65 & 65 & 52 & 117 & 26 & 130 & 0 & 52 & 144 & 26 & 13 & 52 & 0 & 104 & 104 & 117 & 65 & 117 & 104 & 39 & 91 & 130 & 130 & 130 & 0 & 0 \\
156 & 13 & 130 & 39 & 39 & 78 & 91 & 0 & 78 & 79 & 78 & 0 & 0 & 130 & 52 & 156 & 78 & 26 & 104 & 91 & 78 & 13 & 143 & 26 & 130 & 0 \\
104 & 117 & 52 & 117 & 13 & 91 & 130 & 156 & 130 & 117 & 79 & 156 & 0 & 39 & 0 & 91 & 65 & 26 & 26 & 78 & 143 & 143 & 39 & 65 & 65 & 0 \\
78 & 52 & 91 & 52 & 39 & 91 & 0 & 104 & 26 & 143 & 156 & 118 & 0 & 91 & 117 & 78 & 130 & 104 & 143 & 39 & 39 & 78 & 13 & 156 & 104 & 0 \\
0 & 130 & 26 & 117 & 104 & 52 & 39 & 104 & 156 & 104 & 117 & 117 & 1 & 26 & 91 & 0 & 91 & 91 & 91 & 78 & 0 & 156 & 26 & 143 & 39 & 0 \\
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
65 & 26 & 78 & 104 & 52 & 39 & 156 & 78 & 104 & 91 & 52 & 156 & 0 & 78 & 1 & 13 & 26 & 91 & 39 & 13 & 156 & 130 & 13 & 0 & 52 & 0 \\
65 & 117 & 130 & 13 & 65 & 39 & 104 & 91 & 13 & 117 & 91 & 104 & 0 & 104 & 104 & 118 & 13 & 52 & 13 & 52 & 26 & 52 & 13 & 65 & 130 & 0 \\
39 & 130 & 117 & 0 & 156 & 13 & 78 & 104 & 26 & 156 & 39 & 130 & 0 & 52 & 13 & 0 & 1 & 39 & 117 & 117 & 26 & 39 & 65 & 143 & 26 & 0 \\
117 & 52 & 13 & 0 & 130 & 39 & 65 & 143 & 78 & 130 & 117 & 52 & 0 & 156 & 39 & 0 & 0 & 118 & 13 & 13 & 78 & 117 & 26 & 91 & 78 & 0 \\
78 & 130 & 39 & 13 & 117 & 117 & 0 & 13 & 91 & 78 & 39 & 156 & 0 & 78 & 117 & 104 & 0 & 143 & 157 & 26 & 39 & 65 & 52 & 13 & 156 & 0 \\
91 & 130 & 13 & 78 & 39 & 117 & 156 & 143 & 91 & 104 & 143 & 13 & 0 & 156 & 65 & 39 & 91 & 104 & 91 & 157 & 65 & 13 & 13 & 13 & 0 & 0 \\
0 & 104 & 156 & 117 & 104 & 91 & 0 & 104 & 143 & 26 & 13 & 52 & 0 & 65 & 117 & 130 & 39 & 13 & 0 & 13 & 40 & 52 & 91 & 39 & 39 & 0 \\
13 & 13 & 78 & 156 & 78 & 13 & 26 & 117 & 143 & 143 & 143 & 0 & 0 & 104 & 104 & 117 & 52 & 143 & 39 & 0 & 117 & 27 & 143 & 156 & 117 & 0 \\
143 & 91 & 104 & 52 & 130 & 13 & 117 & 52 & 65 & 39 & 130 & 143 & 0 & 13 & 156 & 39 & 130 & 130 & 91 & 78 & 0 & 91 & 92 & 91 & 0 & 0 \\
26 & 0 & 117 & 156 & 130 & 130 & 52 & 39 & 39 & 26 & 156 & 156 & 0 & 65 & 52 & 117 & 52 & 156 & 78 & 39 & 13 & 39 & 52 & 92 & 13 & 0 \\
117 & 78 & 52 & 143 & 13 & 39 & 26 & 26 & 52 & 65 & 104 & 13 & 0 & 91 & 117 & 78 & 117 & 130 & 78 & 0 & 65 & 143 & 26 & 13 & 53 & 0 \\
130 & 117 & 0 & 117 & 117 & 117 & 52 & 0 & 104 & 130 & 39 & 26 & 0 & 0 & 39 & 143 & 52 & 65 & 117 & 130 & 65 & 13 & 65 & 52 & 52 & 1
\end{bmatrix}
$$
As expected, elements of the diagonal are 1 plus a multiple of $p=13$, the other are multiples of $13$. But then $m=pn$ works:
sage: m = p*n
sage: J = matrix.identity(2*p)
sage: M^m == J
True
To get the projective order we are now looking in the tree of divisors of $m$ starting in $m$.
Since i am lazy, i will write the shorter code that produces more work for the machine.
sage: min(d for d in m.divisors() if M^d == (M^d)[0, 0] * J)
878476872
sage: m
878476872
sage: factor(m)
2^3 * 3 * 7 * 13^2 * 30941
Just as a note. The multiplicative order over the field $F$ of the piece $A$ in the block matrix $M$ is:
sage: A.base_extend(F).multiplicative_order()
4826796
sage: _.factor()
2^2 * 3 * 13 * 30941
This is the number that appeared in the OP.
Sorry, there is a typo as the matrix M = block_matrix ([[a, c], [b, d]]). If a is the projective order of M over residual class ring Z/ (169)Z, then, M^{a} = c * mat.identity(26) where c belongs to (Z/(169)Z)^{*}.
But 28392 does not satisfy this definition - see it yourself by running
print( M.change_ring(Zmod(13^2))^28392 ).Actually, It was mentioned in the paper I was reading, so I thought the author was correct. How does one obtain the projective order in this case? I don’t see a way to attach the paper.
Your number 4826796 does not work either.
I know that doesn’t work, but should be the correct projective order?