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!

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}$.
