Why are there $736$ matrices $M\in \mathcal M_2(\mathbb{Z}_{26})$ for which it holds that $M=M^{-1}$? I'm currently trying to introduce myself to cryptography.
I'm reading about the Hill Cipher currently in the book Applied Abstract Algebra. The Hill Cipher uses an invertible matrix $M$ for encryption and the inverse matrix $M^{-1}$ for decryption. The book recommends to use a matrix $M$ for which it holds that $M=M^{-1}$, since we then have the same key for encryption and decryption.
In the book, they use a $2 \times 2$ matrix $M$ over $\mathbb{Z}_{26}$ as example, and state that for there are 736 $2 \times 2$ matrices for which it hold that $M=M^{-1}$.
I'm trying to pick up on as much as possible when reading things, since I find it counter-productive for learning to skip something, when you don't get the theory behind it.
Can someone enlighten to me, as to why it is that there are 736 possible $2 \times 2$ $M=M^{-1}$ matrices and how to find them?
 A: Indeed there turns out to be a more systematic answer. Computing the number $a_n$ of self-inverse matrices over $\mathbb Z_n$ for some $n$ leads to OEIS sequences A066907 and A066947. It turns out that $a_n$ is multiplicative and given by
$$
a_n=c_m\prod_i\left(2+p_i^{2k_i-1}(p_i+1)\right)
$$
for $n=2^m\prod_ip_i^{k_i}$, where the $p_i$ are the odd prime factors of $n$ and
$$
c_m=
\begin{cases}
1&m=0\;,\\
4&m=1\;,\\
28&m=2\;,\\
9\cdot4^{m-1}+32&m\gt2\;.
\end{cases}
$$
In the present case $26=2^1\cdot13^1$, so $a_{26}=4(2+13\cdot14)=736$.
Unfortunately the links in the entry are all broken and the author's email address is obsolete, so the question remains how to derive this.
A: There may be a more systematic way of doing this; here's a semi-systematic way.
Your condition is equivalent to $M^2=I$. For $\displaystyle M=\pmatrix{a&b\\c&d}$ that yields four equations:
$$
\begin{align}
a^2+bc&\equiv1\;,\\
ab+bd&\equiv0\;,\\
ca+dc&\equiv0\;,\\
cb+d^2&\equiv1\;.
\end{align}
$$
The first and last of these have solutions for $a$ and $d$, respectively, if and only if $1-bc$ is a quadratic residue. If so, $a^2\equiv d^2$, and thus $a\equiv\pm d$.
The second and third equation are $(a+d)b\equiv0$ and $(a+d)c\equiv0$. These are fulfilled for $a\equiv-d$, so we have one solution for every pair $b,c$ with $1-bc\equiv0$ or $1-bc\equiv13$ and two solutions for every pair $b,c$ with $1-bc$ a quadratic residue other than $0$ and $13$. I don't know a systematic way of counting these; this code counts $728$.
That leaves the alternative $a\equiv+d$, excluding $a=d=0$ and $a=d=13$ since these have already been counted under $a\equiv-d$. In this case $a+d$ is even and not divisble by $13$, so $(a+d)b\equiv0$ and $(a+d)c\equiv0$ if and only if $b$ and $c$ are both $0$ or $13$. That's $4$ combinations, and for each of them there are two square roots of $1-bc$ that $a=d$ can be, so that makes $8$ additional possibilities, for a total of $728+8=736$.
A: Let $A,B,C$ range over the nonzero elements of ${\mathbb Z}/13$.
Then over ${\mathbb Z}/13$ I get the following self-inverse matrices:
$$\hbox{12 of type} \pmatrix{0&B\cr 1/B&0\cr}$$
$$\hbox{144 of type} \pmatrix{A&B\cr{1-A^2\over B}&-A\cr}$$
$$\hbox{4 of type} \pmatrix{\pm 1&0\cr 0&\pm 1\cr}$$
$$\hbox{12 each of types} \pmatrix{1&0\cr C&-1\cr}, \pmatrix{-1&0\cr C&1\cr},\pmatrix{1&B\cr 0&-1\cr},\pmatrix{-1&B\cr 0&1\cr}$$
That's a total of 208 over ${\mathbb Z}/13$.  If I multiply this by 4 (the number of self-inverse matrices over ${\mathbb Z}/2$), I get 832, not 736.  Where's my mistake?
A: For a $2\times2$ matrix, $M=\begin{pmatrix}a&b\\c&d\end{pmatrix}$, the inverse is given by $M^{-1}=\frac1{ad-bc}\begin{pmatrix}d&-b\\-c&a\end{pmatrix}$. Hence, in order for $M=M^{-1}$, there are two options:
(1) $b=c=0$, then $a=\frac{1}{ad}d=\frac1a$ and $d=\frac1{ad}a=\frac1d$, hence $a^2=d^2=1$, i.e. $a=\pm1$, $b=\pm1$ - so there are $4$ such matrices.
(2) At least one of $b,c$ is non-zero. Then $ad-bc=-1$ and $d=-a$. Two options from here:
(i) If $bc=0$ then $ad=-a^2=-1$, so $a=\pm1$. So there are $25$ options to choose $b\neq0$ and two options to choose $a$. So a total of $50$ matrices with $c=0$. Hence a total of $100$ such matrices with $bc=0$.
(ii) $bc\neq0$. We have $a^2+bc=1$, so $1-bc$ has to be a square $\mod 26$, hence $1-bc\in\{0,1,4,9,16,25,10,23,12,3,22,17,14,13\}$. Since $bc\neq0$, we have $bc\in\{1,23,18,11,2,17,4,15,24,5,10,13,14\}$. There are $26(1-\frac12)(1-\frac1{13})=12$ (Euler's phi) invertible elements in $\mathbb{Z}_{26}$. For any choice of invertible $b$, there are $13$ options for $bc$, each gives one option for $c$, each, in turn gives $2$ options for $a$, except $bc=14$ (since $a^2=13$ implies $a=13$). Hence a total of $12\cdot12\cdot2+12\cdot1\cdot1=300$ such matrices with $b$ invertible.
For $bc\in\{18,2,4,24,10,13,14\}$ there are solutions with $b$ not invertible: each of those has to be checked manually. For example:
$bc=18$ then $b\in\{2,6,9,18,-2=24,-6=20,-9=17,-18=8\}$ (all non-invertible divisors of $18$). Hence there are $8\cdot 2=16$ ($2$ options for $a$) such matrices.
And so on
A: Interesting, this question is very similar (but not completely) to an assignment I have due monday, and which I have tried some of friday and all of saturday to solve (and nearly given up on by now), with no luck. 
I am to answer the exact same question as you, but I am tasked of doing it by "splitting the problem into two cases. First a case where it is modulo 2, then modulo 13". 
I figure that after finding the amount modulo 2 and 13, I should use the Chinese Remainder Theorem to make the result mod 26 again. I cannot figure out though, how I should go about doing it. What I have done so far, is actually something that seems most similar to what joriki suggested (but I have no idea, whether it was the right approach to do in my case). 
What I do to try to solve it, is that i, as instructed, split $det(A) \equiv \pm 1$ mod 26 into: $det(A) \equiv \pm 1$ mod $13$ and $det(A) \equiv \pm 1$ mod $2$. 
Then I have been thinking (and this is the part I think reminded me of joriki's solution) that since $A=A^{-1}$ means that $A^2=I$, I should find all cases where:
$\pmatrix{a&b\\c&d}^2=\pmatrix{1&0\\0&1}$
Which means I have the congruences:
\begin{align}
a^2+bc&\equiv1\\
ab+bd&\equiv0\\
ca+dc&\equiv0\\
cb+d^2&\equiv1\
\end{align}
But that is where the fun ends. I have no clue where to go from this point.... Or if what I have been doing so far, is the correct way of doing it, if I want to show it using the hint given... I hope someone can (and also will help out), and are willing to shed some light on this.
The course in which I have this assignment is an introductory cryptology course, so I am not familar (yet, I find math compelling, especially in a cryptology setting) with Euler's phi, OEIS sequences and quadratic residue terms, so I haven't been able to turn them around and give me enough inspiration to solve my own subproblem here.
Sorry if I am hijacking your question. I made a separate question at first, but it was considered a duplicate to this. And I agree with that. This is a subquestion about how to solve your question in a specific way. I hope the answer (I hope I get one) might be interesting to you too, since you seem to be introducing yourself to cryptography as well. :)
