Symmetric matrix factorization of a $2 \times 2$ symmetric matrices Let S be a symmetric matrix such that
$$S= \begin{bmatrix}
A & B\\
B & C
\end{bmatrix}, \text{where A, B, C $\in$ $\mathbb{Z_P},$}$$
 where p is prime.
I am tasked to get the solution(s) of the symmetric matrix factorization of $2 \times 2$ symmetric matrices of the forms,
$$S = X^2(\mod p)$$ 
and 
$$S = XY(\mod p),$$
where $X$ $\neq$ $Y$ but are both symmetric.
For example, under $\mathbb{Z_3}$,
$$S = \begin{bmatrix}
2 & 2\\
2 & 1
\end{bmatrix} = \begin{bmatrix}
1 & 2\\
2 & 0
\end{bmatrix} \begin{bmatrix}
1 & 2\\
2 & 0
\end{bmatrix}(\mod 3) = X^2$$
and
$$S = \begin{bmatrix}
1 & 2\\
2 & 2
\end{bmatrix} = \begin{bmatrix}
2 & 0\\
0 & 1
\end{bmatrix} \begin{bmatrix}
2 & 1\\
2 & 2
\end{bmatrix}(\mod 3) = XY$$
I know that these pairs are not the only possible factors of the given matrix S. I've also done programming the matrix multiplication which is the inverse of matrix factorization because I have no idea how to start with matrix factorization.
Any suggestions will be much appreciated. Thanks in advance!
 A: $S$ is symmetric and hence disagonalisable. First, we find the eigenvalues and eigenvectors of $S$. The characteristic polynomial is $\chi_S(t)=|S-tI|=(A-t)(C-t)-B^2$, which has roots $\lambda_1=\frac{1}{2}\left(\sqrt{A^2-2AC+4B^2+C^2}+A+C\right), \lambda_2=\frac{1}{2}\left(-\sqrt{A^2-2AC+4B^2+C^2}+A+C\right)$. Let $v=(v_1, v_2)^t, u=(u_1, u_2)^t$ be their corresponding eigenvalues. Then we have $Sv=\lambda_1 v, Su=\lambda_2 u$. A tedious calculation gives a possible solution $v=\left(-\frac{-A+C+\sqrt{A^2-2AC+4B^2+C^2}}{2B}, 1\right)$ and $u=\left(\frac{A-C+\sqrt{A^2-2AC+4B^2+C^2}}{2B}, 1\right)$. Let $P$ be the matrix with columns $v, u$ and $D=\operatorname{diag}(\lambda_1, \lambda_2)$. Then $S=PDP^{-1}$. Hence, letting $D'=\operatorname{diag}\left(\sqrt{\lambda_1}, \sqrt{\lambda_2}\right)$ we get that $D'^2=D$ and hence $S=PD'^2P^{-1}=(PD'P^{-1})(PD'P^{-1})$. Now setting $X=PD'P^{-1}$ we get $S=X^2$. Of course, we can set $X=PD'$ and $Y=P^{-1}$ to get $S=XY$. Note I assumed $A^2-2AC+4B^2+C^2, \lambda_1, \lambda_2$ turn out to be squares in $\mathbb{Z}_p$, and I'm uncertain on how to prove this; I'm not even sure if it is true. Suggestions in the comments would be appreciated.
A: First, if $S$ is diagonal then this is fairly straightforward. We distinguish three cases:


*

*If $S=0$ then there are $(p-1)^2(p+1)$ pairs $(X,Y)$ of nonzero matrices such that $S\equiv XY\pmod{p}$, and $2p^4-1$ such pairs $(X,Y)$ with either $X=0$ or $Y=0$.

*If $B=0$ and $A=C\neq0$ then $S=AI$, and hence the solutions to
$$S\equiv XY\pmod{p},$$
are precisely the pairs $(X,Y)$ where $X$ is invertible and $Y=AX^{-1}$. In particular there are $(p-1)^2p(p+1)$ solutions $(X,Y)$.

*If $B=0$ and $A\neq C$ then both $X$ and $Y$ must be diagonal and so
$$A\equiv x_{11}y_{11}\pmod{p}\qquad\text{ and }\qquad B\equiv x_{22}y_{22}\pmod{p}.$$
There are $(p-1)^2$ solutions if $A\neq0$ and $C\neq0$, and $(p-1)(2p-1)$ solutions if $AC=0$.
From here one we assume that $B\neq0$. If $A=C=0$ then for the matrix $R:=\binom{0\ 1}{1\ 0}$ we have
$$SR=\begin{pmatrix}
0&B\\B&0
\end{pmatrix}
\begin{pmatrix}
0&1\\1&0
\end{pmatrix}
=\begin{pmatrix}
B&0\\0&B
\end{pmatrix},$$
so $SR$ is diagonal. The above tells us the solutions $(X,Y)$ for $SR$, and hence because $R^2=I$ we have
$$S\equiv XYR\pmod{p},$$
so the solutions for $S$ are the pairs $(X,YR)$ where $SR\equiv XY\pmod{p}$.
If $A$ and $C$ are not both zero then wlog $A\neq0$. For the matrix $T:=\binom{1\ \hphantom{-}B}{0\ -A}$ we have
$$T_A^{\intercal}ST_A=\begin{pmatrix}
1&0\\B&-A
\end{pmatrix}\begin{pmatrix}
A&B\\B&C
\end{pmatrix}
\begin{pmatrix}
1&B\\0&-A
\end{pmatrix}
=\begin{pmatrix}
A&0\\0&A|S|
\end{pmatrix},$$
where $|S|=\det S$, so $T_A^{\intercal}ST_A$ is diagonal. Again the argument above tells us the solutions for $T_A^{\intercal}ST_A$. Note that $T_A$ is invertible because $A\neq0$, and so
$$S\equiv T_A^{-\intercal}XYT_A^{-1}\pmod{p},$$
and hence the solutions for $S$ are the pairs $(T_A^{-\intercal}X,YT_A^{-1})$ where $T_A^{\intercal}ST_A\equiv XY\pmod{p}$.
