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$?
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.