Square root of Matrix over $F_p$ Let $M \in F_p^{n\times n}$ be a square matrix of dimension $n\times n$ with entries over $F_p$.
Then we want to obtain $M$ from $M^2=M\cdot M$. Initially I think of the old diagonalisation procedure of solving the square root of a matrix by:
$$M^2 = P\cdot D \cdot P^{-1}$$
$$M = P\cdot D^{\frac{1}{2}} \cdot P^{-1}$$
Where $D$ is a diagonal matrix with the eigenvalues as entries defined over $F_p$. So $2^{-1} \pmod p$ must be a unit on $F_p^{*}$, this is $2$ must be an unit on $F_p^{*}$ so it has an inverse.
Moreover, when taking $F_2$ as the field, the diagonalisation method doesn't yield the square root matrix as $2$ is not an unit on $F_2^{*}$.
I know that if $M^k \in GL(n,F_2)$ then $\gcd(k,\vert GL(n,F_2) \vert)$ must be equal to $1$ so $M^{k\cdot d}=M$ where $d$ is the multiplicative inverse of $k$ modulo the order of the given general linear group. It works also with the multiplicative order of $M$ in $GL$.
But what if $k=2$ and the field is $F_2$? Is there any method for computing the square root matrix of $M^2$ over $F_2$?
 A: Depends.
You already start with a rather strong assumption, that your matrix $M^2$ can be diagonalized. That is wrong in general and only true for some matrices.
But let's assume that you can write
$$M^2 = PDP^{-1}.$$
Then every diagonal entry of $D$ is either $0$ or $1$, and both of these have unique square roots. In fact, you get $D^{1/2} = D$ and thus $M = M^2$.
That is, of course, a boring case, but once again shows that your assumption that $M$ is diagonizable is already very restrictive; in the case of the binary field, it means that $M^k = M$ for all powers $k$.
In general it seems you are slightly confused with some things. To be able to get a unique $k$-th root, you need that the map
$$f : \mathbb{F}_p \to \mathbb{F}_p, x \mapsto x^k$$
is a bijection. This is equivalent to saying that the equation $x^k = 1$ has a unique solution in $\mathbb{F}_p$. According to this answer here, that is equivalent to $gcd(k,p-1) = 1$. That means that you can take unique square roots only over $\mathbb{F}_2$, over all other finite prime fields there are elements that do not have a square root.
So it turns out that the case of $\mathbb{F}_2$ is the easy one here, if you assume that such a matrix $D$ exists.
For very general $M$ it seems finding a root is NP complete, so I don't think you will find an efficient way to do it.
A: As Dirk wrote , your post is very unclear. 
We study the following
$\textbf{Problem.}$ Let $A\in M_n(F_2)$; what are the solutions $R\in M_n(L)$ (if they exist) of the equation $R^2=A$ ? (here $L$ is an algebraic extension of $F_2$.
We consider the (unique) Chevalley-Jordan decompositions of $A,R$: $A=D+N,R=\Delta +M$ where $D,\Delta$ are diagonalizable on an extension of $F_2$ and $N,M$ are nilpotent and $DN=ND,\Delta M=M\Delta$.
$\textbf{Remark}$. Note that $D,N\in M_n(F_2)$ can be effectively calculated for every matrix $A$. The unknowns are $\Delta,M$.
Thus $R^2=\Delta^2+M^2=A=D+N$, that implies $(1)$ $\Delta^2=D$ and $(2)$ $M^2=N$.


*

*We give two examples for $(1)$ (there exist always solutions in $\Delta$)


i) Assume that $\chi_D(x)=p(x)=x^5+x^2+1$, a polynomial that is irreducible over $F_2$. Their roots are $a,a^2,a^4,a^8,a^{16}$ where $a\in K=F_2[x]/p(x)$, an algebraic extension of $F_2$ of degree $5$ ($a=x$ in the quotient).
Then $D=Pdiag(a,a^2,a^4,a^8,a^{16})P^{-1}$ (where $P\in GL_n(K)$) and the UNIQUE solution in $\Delta$ is $\Delta=Pdiag(a^{16},a,a^2,a^4,a^8)P^{-1}$.
EDIT. Note that $D,\Delta$ are similar. In fact, $\Delta\in M_5(F_2)$. Indeed, $D$ is similar in $M_5(F_2)$ to the companion matrix $D'=\begin{pmatrix}0&0&0&0&1\\1&0&0&0&0\\0&1&0&0&1\\0&0&1&0&0\\0&0&0&1&0\end{pmatrix}$ and its square root is $\Delta'=\begin{pmatrix}1&1&1&0&0\\1&1&1&1&0\\0&0&0&1&1\\1&0&0&0&1\\1&1&0&0&0\end{pmatrix}$.
More precisely $\Delta=D^{16}$ (cf. the associated eigenvalues).
ii) Note that, if $D=0$, then $\Delta=0$.


*

*and two examples for $(2)$ (There don't exist always solutions in $M$).


i) No solutions when $N=J_n$ (the Jordan block of dimension $n$).
ii) Let $N=diag(J_2,J_2)$. Then there are $12$ solutions in $M_4(F_2)$, including the following $4$:
$M=\begin{pmatrix}0&a&1&b\\0&0&0&1\\0&1&0&a\\0&0&0&0\end{pmatrix}$.
