Let $p$ be an odd prime and $\mathbb Z_p$ be the prime field of order $p$. Consider the matrix ring $R=M_n(\mathbb Z_p)$. Is there any method to count the solutions of the equation (in the ring $R$)
$$X^2=I.$$ Where $I$ is the identity matrix?
Let $p$ be an odd prime and $\mathbb Z_p$ be the prime field of order $p$. Consider the matrix ring $R=M_n(\mathbb Z_p)$. Is there any method to count the solutions of the equation (in the ring $R$) $$X^2=I.$$ Where $I$ is the identity matrix? 


Denote $k:=Z_p$. If $p\ne2$, the matrix $X$ is in onetoone correspondence with a decomposition $k^n=E_+ \oplus E_$, where $E_\pm$ is the eigenspace associated with the eigenvalue $\pm1$. Given the dimension $m$ of $E_+$ ($nm$ for $E_$), these decompositions are in onetoone correspondence with the bases of $k^n$, quotiented by the action of $GL_m\times GL_{nm}$. Therefore, their number is $$\frac{GL_n(Z_p)}{GL_m(Z_p)\cdotGL_{nm}(Z_p)}=p^{m(nm)}\frac{(p1)\cdots(p^n1)}{(p1)\cdots(p^m1)(p1)\cdots(p^{nm}1)}.$$ Summing up over all the possible $m$'s, the number of solutions to $X^2=I$ is $$\sum_{m=0}^np^{m(nm)}\frac{(p^{m+1}1)\cdots(p^n1)}{(p1)\cdots(p^{nm}1)}.$$ 


Because $p$ is different from $2$, the solutions are the same as direct sum decompositions $\mathbb{Z}_p^{\oplus n} = E_{+1} \oplus E_{1}$ as $\mathbb{Z}_p$vector spaces. You can index these by the dimensions, say $r$ and $nr$. For each, the number of solutions is $$ \frac{(p^n1)(p^np) \cdot \dots \cdot (p^np^{r1})}{(p^r1)(p^rp) \cdot \dots \cdot (p^rp^{r1})} \cdot p^{r(nr)} $$ So the final answer is the sum over $r$ from $0$ to $n$ of this factor. 


There are $\displaystyle \sum_{i=0}^n \frac{N(n)}{N(i)N(ni)}$ solutions to your problem, where $N(r)$ is the number of elements in $GL_r(\mathbb{F}_p)$. Any such $X$ is invertible, and $X\in GL_n(\mathbb{F}_p)$ is a solution if and only if it is conjugate in $GL_n(\mathbb{F}_p)$ to a diagonal matrix with the first $i$ entries equal to $1$ and last $ni$ entries equal to $1$. This explains the number of summands. Now $g\in GL_n(\mathbb{F}_p)$ commutes with such diagonal matrix if and only if it preserves both eigenspaces. Hence, there are $\frac{N(n)}{N(i)N(ni)}$ elements in the conjugacy class. The formula can be made even more explicit by substituting a formula for $N(r)$. 

