I just analysed this equation for real matrices $$ A^2=\begin{pmatrix}a&b\\c&d\end{pmatrix}^2=I $$ From the main diagonal of $A^2$ we must have $a^2+bc=bc+d^2=1$ showing that $d=\pm a$.
CASE 1: If $b=0$ we have $a^2=d^2=1$ so then $c=0$. If $c=0$ the same argument renders $b=0$.
CASE 2: Otherwise assume $b,c\neq 0$. Then if $d=a$ the opposite diagonal of $A^2$ shows that $2ab=2ac=0$. Thus $a=d=0$ and the equation $a^2+bc=1$ yields $bc=1$.
CASE 3: Now, if $d=-a$ we have $ab-ab=ac-ac=0$ which is always true, so then only the first equation $a^2+bc=1$ matters.
So we have three sets of solutions:
- CASE 1: gives the four solutions $b=c=0$ and $a,d\in\{-1,1\}$.
- CASE 2: gives an infinite family of solutions $a=d=0$ and $b=x,c=x^{-1}$ with $x\in\mathbb R\setminus\{0\}$.
- CASE 3: gives a last family of solutions: choose $(x,y)\in\mathbb R\times\left(\mathbb R\setminus\{0\}\right)$ and define $a=x,d=-x$ and $b=y,c=(1-x^2)y^{-1}$.
I find it interesting how the solutions divide nicely into advancing stages of first discrete, then one-dimensional, and finally two-dimensional parametrisations in this analysis. I wonder if that is purely a coincidence ...
QUESTION: What would a similar analysis of $A^2=I$ where we are dealing with $3\times 3$ matrices look like? The case of $3\times 3$ looks rather complicated already, so I would hardly dream of considering the general case of $n\times n$ matrices.