Number of self-inverse matrices over prime field Regarding the cryptosystem known as the Hill Cipher, my textbook by Douglas R. Stinson has an exercise asking you to find the number of involutory keys for $m=2$ over $\mathbb Z_{26}$. This means that we are to find the number of $2\times 2$ matrices over $\mathbb Z_{26}$ that satisfies $K=K^{-1}$. This can be considered a basic mathematical problem regardless of the details of the Hill Cipher cryptosystem.

I have found the number to be $736$ by applying the easily proven result
$$
K=K^{-1}\implies\det(K)=\pm1\pmod{26}
$$
stated in the previous exercise, and the general formula
$$
\begin{pmatrix}
a&b\\c&d
\end{pmatrix}^{-1}
=\det(K)^{-1}
\begin{pmatrix}
d&-b\\-c&a
\end{pmatrix}
$$
stated elsewhere in the textbook and hinted at in the exercise in question.

But here is my problem:

This approach seemed far from general, dividing into cases and subcases regarding whether $\det(K)=1$ or $\det(K)=-1$ and regarding the profiles of $a,b,c$ modulo $2,13$ respectively (Chineese Remainder Theorem [CRT]). This lead to conclusions equivalent to saying that for $\det(K)=-1$ the equation $K=K^{-1}$ has $4$ solutions over $\mathbb Z_2$ and $182$ solutions over $\mathbb Z_{13}$. By [CRT] this provides $4\cdot182=728$ solutions over $\mathbb Z_{26}$ for the subcase $\det(K)=-1$. The case $\det(K)=1$ provided $8$ matrices in my analysis.

So my question is:

Is there a general approach for finding the number of solutions to the matrix equation $K=K^{-1}$ over a prime field $\mathbb Z_p$ for matrices of arbitrary size $m\times m$?

 A: This is closely related to the number of invertible matrices over a finite field.
If $M$ is an $m \times m$ matrix over $\mathbb{F}_p$ then $M^2 - I = 0$ implies that the minimal polynomial of $M$ divides $t^2 - 1$.
Then if this is the case every eigenvalue of $M$ is either $1$ or $-1$. So we can put $M$ into Jordan normal form 
$$
M = P^{-1} J P,
$$
where $P$ is invertible and $J$ is in Jordan normal form with only $1$ and $-1$ along the diagonal.
Note though that if $M^2 = I$ then 
$$
P^{-1}JP P^{-1}JP = I, \text{ so } J^2 = I.
$$
It is easy to see that over any field of positive characteristic that this implies that all Jordan blocks are of size 1.
So then the number of diagonal matrices with $1$ or $-1$ along the diagonal is $2^m$, but we want to exclude the case when all the entries are $1$ or all entries are $-1$ because then $P^{-1} J P = \pm I$, and over a field of characteristic 2 this is the only possibility so we will have to deal with that separately, so we have $2^m - 2 $ possible diagonal matrices.
We then have to multiply this number by the number of invertible $m \times m$ matrices over $\mathbb{F}_p$ and it is pretty well known that this number is
$$
\prod_{i=0}^{m-1}(p^m - p^i).
$$
See https://www.math.wisc.edu/~ddrake/pdf/sl_n_q.pdf if you need details.
So this gives us a grand total of
$$
\bigg{[}(2^m -2)\prod_{i=0}^{m-1}(p^m - p^i)\bigg{]} + 2
$$
matrices that are their own inverse over $\mathbb{F}_p$. The plus two at the end is the identity matrix and its negative, which you might want to discount (you probably wouldn't want to discount its negative though). Over characteristic $2$ this is not valid as $1 = -1$ so the only diagonal matrix with $1$ or $-1$ along the diagonal is the identity matrix.
Now to show these matrices are all distinct we suppose that $P^{-1}DP = Q^{-1}D'Q$ where $D, D'$ are diagonal and $P, Q$ are invertible. If these matrices are equal then in particular they have the same eigenvalues and in the same order in the Jordan normal form, that is $D = D'$. If you want more convincing of this, note that these matrices are equal, so have the same eigenvectors, with the same eigenvalues, so $D =D'$.
Also their eigenvectors are the same, so as the columns of $P$ are a basis of  eigenvectors, these must be the same in $P$ and $Q$ UP TO MULTIPLICATION BY A NON-ZERO SCALAR. So in our initial sum we have counted too many, we need to take multiplication by a scalar into account. But this just means we need to divide by $(p-1)^m$ giving us a grand total of
$$
\frac{(2^m - 2)\prod_{i=0}^{m-1}(p^m - p^i)}{(p-1)^m}
$$
plus the identity matrix and its negative.
Note that if we impose the condition det$(M) = -1$ then this gives us half of this number, over odd characteristic. If $\lambda_1, \ldots \lambda_m$ are the eigenvalues of $M$ then 
$$
\text{det}(M) = \prod_{i=1}^m \lambda_i.
$$
so the condition det$(m) = -1$ fixes the last eigenvalue, so there are only $2^{m-1} - 1$ possible diagonal matrices, giving us (when $p =13$) $182$ possible matrices which are self inverse. Note that as $m$ is even we don't include the negative of the identity matrix.  
