how many $2\times2$ matrices exist with this condition? Let $p=10007$ which is prime. I want to find the number of matrices X of $2\times2$ dimension with elements from $\mathbb{Z}_p$ for which $X^2\equiv I$ (mod $p$). Where $I$ - identity matrix. How can I solve this problem? I tried to write down the square of matrix so that I get the system of congruences, but I can't figure out what to do after that, and I don't think this is the right way to solve this.
 A: Let's say $X=\begin{bmatrix} a&b\\c&d\end{bmatrix}$, after computing $X^2$, we get the following system of congruences
$$\begin{cases} a^2+bc\equiv 1\mod p\\b(a+d)\equiv 0\mod p\\c(a+d)\equiv 0\mod p\\d^2+bc\equiv 1\mod p\end{cases}$$
The $2$nd and $3$rd equations are simple to solve, so we will do casework on the those. We either have $b,c\equiv 0\mod p$ or $a+d\equiv 0\mod p$. There is also some overlap between these two cases, so we have to watch out for overcounting
Case 1: $b,c\equiv 0\mod p$
Then we have
$$\begin{cases} a^2\equiv 1\mod p\\d^2\equiv 1\mod p\end{cases}$$
There are clearly $2$ solutions for each of the congruences ($a,d\equiv \pm 1\mod p$), for a total of $4$ solutions. Since the $2$ solutions when $a\equiv -d\mod p$ overlap with the second case, we will exclude those for a total of $\boxed{2}$.
Case 2: $a+d\equiv 0\mod p$
Then we have
$$bc\equiv 1-a^2\equiv 1-d^2\mod p$$
We will split this case into sub cases depending on the values of $a,d$.
Case 2a: $a^2,d^2\equiv 1\mod p$
The solutions for $a,d$ are $a\equiv \pm 1\mod p,d\equiv \mp 1\mod p$. Hence there are $2$ solutions for $a,d$.
To find the number of solutions for $b,c$, we need to find the number of solutions to
$$bc\equiv 0\mod p$$
This only happens if at least $1$ of $b,c$ is equivalent to $0\mod p$. We must be careful to not double count the solutions of $b,c\equiv 0\mod p$. We get a total of $p+p-1=2p-1$ solutions for $b,c$.
Hence, there are $2(2p-1)=\boxed{4p-2}$ solutions for this subcase.
Case 2b: $a^2,d^2\equiv k\not\equiv 1\mod p$
Since $a,d\not\equiv -1,1\mod p$, there are $p-2$ ordered pairs to choose for $(a,d)$. Each of these generate a $k$ that satisfies $1-k\not\equiv 0\mod p$.
To find the solutions for $b,c$, these satisfy
$$bc\equiv 1-k\mod p$$
Since $1-k\not\equiv 0\mod p$, we have that
$$gcd(b,p)=1\implies c\equiv (1-k)b^{-1}\mod p$$
Hence, there is a unique solution for all $b\in[1,p-1]$.
The total in this subcase is $(p-1)(p-2)=\boxed{p^2-3p+2}$

Totaling all of these up, the number of solutions for $X$ is
$$2+(4p-2)+(p^2-3p+2)$$
$$\boxed{\boxed{p^2+p+2}}$$
This formula should only work when $p>2$. This is because throughout the process we assumed $1\not\equiv -1\mod p$. Upon substituting $p=10007$, we get the total number of solutions is $\boxed{100150058}$.
A: As in @ThomasAndrews' comment, there is indeed a systematic (not-so-computational) way to approach this, maybe less lending itself to computational errors and such.
So, for a two-by-two matrix with $x^2=1$, certainly $x^2-1=0$, so either $x-1=0$, or $x+1=0$, or neither holds, but $x^2-1=0$.
In the first two cases, $x$ is $\pm$ the identity matrix, that is, has either eigenvalue $1$ (with multiplicity $2$) or eigenvalue $-1$ (with multiplicity $2$).
In the third case, necessarily $x$ has eigenvalues $1$ and $-1$. By standard linear algebra (!!!), this means that there is an invertible two-by-two matrix $A$ such that
$$
x \;=\; Ax_oA^{-1}
$$
with
$$
x_o \;=\; \pmatrix{1 & 0 \cr 0 & -1}
$$
Two matrices $A$ and $B$ produce the same $x$ if and only if
$$
x_o \;=\; (B^{-1}A)x_o(B^{-1}A)^{-1}
$$
That is, $B^{-1}A$ is in the normalizer of $x_o$. We can check that this normalizer is a subgroup of the two-by-two invertible matrices. For $g=\pmatrix{a & b\cr c & d}$, the condition $gx_og^{-1}=x_o$ (multiplied out) is equivalent to $b=c=0$, if $-1\not=+1$, which does hold for $p\not=2$. That is, for $p\not=2$, the normalizer consists of diagonal matrices.
Thus, the collection of matrices $x$ with eigenvalues $1$ and $-1$ is in bijection with the coset space of two-by-two invertible matrices by the subgroup of diagonal invertible matrices.
The number of diagonal matrices is $(q-1)^2$, while the cardinality of two-by-two invertible matrices is $(q^2-1)(q^2-q)$. Thus, the number of such matrices $x$ is the quotient, namely, $(q+1)q$.
Thus, adding the identity and negative identity, we get $(q+1)q+1+1=q^2+q+2$. (Matching Alan Abrahams more concrete computation.)
