Taking Square Roots of Matrices over $\mathbb{Z}/n\mathbb{Z}$ Is it easy (computationally) to take square roots of matrices over $\mathbb{Z}/n\mathbb{Z}$, if you know the factorization of $n$?
If the matrix is diagonalizable, then does diagonalizing and taking the square roots work?  What if the matrix isn't diagonalizable, or if the diagonal has entries which aren't squares?  What does this imply about square roots of the matrix?
Sorry for not phrasing this question well originally.  I am interested in finding ANY square root.  But finding a random square root is even better.  Also, you can imagine the input to be the square of a random matrix.  
 A: We consider the following problem. 
Let $A\in M_k(\mathbb{Z}_n)$ be considered as unknown and $B=A^2$ is given.
Can we obtain a matrix $X\in M_k(\mathbb{Z}_n)$ s.t. $(*)$ $X^2=B$ ?
Assume that the decomposition of $n$ in primes is squarefree: $n=p_1\cdots p_t$. Then it suffices to solve the equations $(*)$ successively in $\mathbb{Z}_{p_1},\cdots,\mathbb{Z}_{p_n}$. After, we lift the results in $\mathbb{Z}_n$. In the sequel, $n=p$ is a prime.
$\textbf{Proposition}$. If $p>2$, $A$ has no non-zero opposite eigenvalues and $dim(\ker(B))\leq 1$, then $A\in\mathbb{Z}_p[B]$.
$\textbf{Proof}$. Up to a change of basis in $K$ (an algebraic extension of $\mathbb{Z}_p$), we may assume that $A=diag(\lambda_1I_{i_1}+N_1,\cdots)$ where the $(\lambda_i)$ are distinct and the $(N_i)$ are nilpotent and if $\lambda_j=0$, then $i_j=1$; consequently, the eigenvalue $0$ is not a problem and, in the sequel, we assume that the $(\lambda_i)$ are $\not= 0$. Note that $B=diag(\lambda_1^2I_{i_1}+2\lambda_1N_1+N_1^2,\cdots)$ where the $\lambda_i^2$ are distinct.
Since $p>2$, $B_1=\lambda_1^2I_{i_1}+2\lambda_1N_1+N_1^2$ is similar to $\lambda_1^2I_{i_1}+N_1$ and the matrices that commute with $B_1$ are polynomials in $B_1$; then $A_1=\lambda_1I_{i_1}+N_1$ is a polynomial in $B_1$, and consequently, $A=r(B)$ where $r\in K[x]$ has  a degree $<d$, the degree of the minimal polynomial of $B$.
${I_k,B,\cdots,B^{d-1}}\subset M_k(\mathbb{Z}_p)$ is a free system over $K$; since $A\in M_k(\mathbb{Z}_p)$, the coefficients of $r$ are in $\mathbb{Z}_p$ and we are done.  $\square$ 
Conclusion. Under the above hypothesis (satisfied when we randomly choose $A$), we no longer have $k^2$ unknowns but only $k$. 
If you don't want to consider finite fields extensions by passing to the quotient, you can do the following simple reasoning when $n$ is small
EXAMPLE. $n=pq,p=47,q=53,k=4$.     
 
The unknowns are $a,b,c,d$ and we calculate $C=(aI_4+bB+cB^2+dB^3)^2-B$.
1) We solve $C=0$ in $M_4(\mathbb{Z}_p)$ by varying $a,b,c,d$ in $[[1..p]]$; we stop when we obtain the first solution.

2) The same method in $M_k(\mathbb{Z}_q)$ gives 

3) The final lift gives a solution $\not= A$.

Time of calculation 2'40"
