$P(AB=BA)$ , $A,B\in M_{3x3}(\mathbb Z/p\mathbb Z)$ Let $A,B\in M_{3x3}(\mathbb Z/p\mathbb Z)$ ($p$ a prime number). Find the probability $P$ that $AB=BA$ that is $P(AB=BA)$
$$A=\begin{pmatrix} a_{11} & a_{12} & a_{13} \\ a_{21} & a_{22} & a_{23} \\ a_{31} & a_{32} & a_{33} \\ \end{pmatrix} $$
$$B=\begin{pmatrix} b_{11} & b_{12} & b_{13} \\ b_{21} & b_{22} & b_{23} \\ b_{31} & b_{32} & b_{33} \\ \end{pmatrix} $$
Please I would really appreciate if you can help me with this problem. Any ideas or suggestions would be highly appreciated :)
 A: Let $P_p$ be the required probability. The couples $\{(A,uI_3+vA+wA^2)|A\in M_3\}$ are some solutions ; moreover, "generically", they are THE solutions. Then $P_p\approx \dfrac{p^9p^3}{p^{18}}=\dfrac{1}{p^6}$ ; moreover $P_p\geq 1/p^6$. Numerical experiments (for $p=2,3,5,7$) seem to "show" that $1/p^6\leq P_p\leq 2/p^6$.
A: This is perhaps not very satisfying, but brute-forcing the computation for $p = 2$ in Maple with the below code gives that there are $7456$ commuting (ordered) pairs of matrices in $M_{3 \times 3}(\mathbb{Z} / 2\mathbb{Z})$ (out of $(2^9)^2$ total pairs), so the probability in that case is
$$\frac{7456}{(2^9)^2} = \frac{233}{8192} = 0.028442\ldots,$$
which agrees with loup blanc's experiments.
p := 2;
M := [seq(seq(seq(seq(seq(seq(seq(seq(seq(Mod(p, Matrix(3, (i, j) -> m||i||j), integer), m33 = 0..(p-1)), m32 = 0..(p-1)), m31 = 0..(p-1)), m23 = 0..(p-1)), m22 = 0..(p-1)), m21 = 0..(p-1)), m13 = 0..(p-1)), m12 = 0..(p-1)), m11 = 0..(p-1))]:
k := 0;
for i from 1 to nops(M) do
    for j from 1 to nops(M) do
        if op(convert(Mod(p, M[i].M[j] - M[j].M[i], integer[]), set)) = 0 then k := k + 1 end if;
    end do;
end do;
k;
k/

(Incidentally, if someone knows a better way to generate a list of all matrices in $M_{3 \times 3}(\mathbb{Z} / p \mathbb{Z})$ I'd be grateful to learn how.)
A: Edit on 10/13/2016
The number of similarity classes $|I|$ of $\mathrm{GL}_3(\mathbb{F}_p)$ below easily follows from this answer. In fact, without having to consider $8$ different types, we have
$$
|I|=\sum_{\lambda \in P_3} p^{\ell(\lambda)}
$$
where $P_3$ is the set of partition of $3$, and $\ell(\lambda)$ is the number of parts of $\lambda$. 
Then 
$$
|I|=p^3+p^2+p.
$$
An Upper bound for  $p\geq 5$
Let $A\in M_3(\mathbb{Z}/p\mathbb{Z})$. Denote by $\mathbb{F}_p=\mathbb{Z}/p\mathbb{Z} $ and $\mathbb{F}_p[x]$ the polynomial ring over $\mathbb{F}_p$. We write $M^A=\mathbb{F}_p^3$ the finitely generated module over $\mathbb{F}_p[x]$ where $x \cdot v = A \cdot v$ for any $v\in\mathbb{F}_p^3$. It is well-known that the dimension of the space of commuting matrices $C_A = \{ B \in M_3(\mathbb{Z}/p\mathbb{Z}) | AB=BA\}$ depends only on the similarity class of $A$. The similarity class of $A$ falls into one of the following eight types:


*

*$aaa$-type: $M^A \simeq \frac{\mathbb{F}_p[x]}{(x-a)}\oplus \frac{\mathbb{F}_p[x]}{(x-a)}\oplus \frac{\mathbb{F}_p[x]}{(x-a)}$ where $a\in\mathbb{F}_p$.

*$aa^2$-type: $M^A \simeq \frac{\mathbb{F}_p[x]}{(x-a)}\oplus \frac{\mathbb{F}_p[x]}{(x-a)^2}$ where $a\in\mathbb{F}_p$.

*$a^3$-type: $M^A \simeq \frac{\mathbb{F}_p[x]}{(x-a)^3}$ where $a\in\mathbb{F}_p$.

*$C$-type: $M^A \simeq \frac{\mathbb{F}_p[x]}{(C(x))}$ where $C(x)\in\mathbb{F}_p[x]$ is a monic cubic irreducible polynomial over $\mathbb{F}_p$. 

*$abb$-type: $M^A \simeq \frac{\mathbb{F}_p[x]}{(x-a)}\oplus \frac{\mathbb{F}_p[x]}{(x-b)}\oplus \frac{\mathbb{F}_p[x]}{(x-b)}$ where $a, b\in\mathbb{F}_p$ and $a\neq b$.

*$ab^2$-type: $M^A \simeq \frac{\mathbb{F}_p[x]}{(x-a)}\oplus \frac{\mathbb{F}_p[x]}{(x-b)^2}$ where $a, b\in\mathbb{F}_p$ and $a\neq b$.

*$aB$-type: $M^A \simeq \frac{\mathbb{F}_p[x]}{(x-a)}\oplus \frac{\mathbb{F}_p[x]}{(B(x))}$ where $a \in\mathbb{F}_p$ and $B(x)$ is a monic quadratic irreducible polynomial over $\mathbb{F}_p$.  

*$abc$-type: $M^A \simeq \frac{\mathbb{F}_p[x]}{(x-a)}\oplus \frac{\mathbb{F}_p[x]}{(x-b)}\oplus \frac{\mathbb{F}_p[x]}{(x-c)}$ where $a,b,c\in\mathbb{F}_p$ are distinct.
Let $\{I_1,I_2,\ldots, I_8\}$ be the set of the above eight types similarity classes, and let $I=\cup_{j=1}^8 I_j$ be the set of all distinct similarity classes.  Define $G=GL_3(\mathbb{F}_p)$, and $G_A=\{H\in G| AH=HA\}$. By the orbit-stabilizer formula, 
$$
|\{(A,B)\in M_3(\mathbb{F}_p)^2 | AB=BA\}|=\sum_{A\in I} |G/G_A| |C_A|\leq \sum_{A\in I} |G| \frac{|C_A|}{|C_A|-3p^{\dim C_A -1}}.
$$
For the last inequality, let $A\in I_j$ and $\dim C_A = m$. Let $\{I=B_1, \ldots, B_m\}$ be a basis for $C_A$. Consider 
$$
\det (c_1 B_1 + \cdots + c_m B_m) = 0.  \ \ \ (*)
$$
For any given $c_2,\ldots, c_m$, the equation $(*)$ is a cubic equation in $c_1$. Therefore, the number of solution to $(*)$ is at most $3p^{m-1}$. 
We remark that with a lot of effort, we will be able to find explicit formula for $|G_A|$. This is certainly more difficult than finding $|C_A|$.
Note that 
$$
1\leq\frac{|C_A|}{|G_A|}\leq\frac{|C_A|}{|C_A|-3p^{\dim C_A -1}}=\frac1{1-\frac3p}=\frac p{p-3}.
$$
The number of elements in eacy $I_j$ can be easily found:
$$|I_1|= |I_2|=  |I_3|=p, \ |I_4|= \frac{p^3-p}3, \ |I_5| =|I_6|=p(p-1), \ |I_7|=\frac{p(p^2-p)}2,$$
$$|I_8|=\frac{p(p-1)(p-2)}6.$$
Therefore, 
$$
|\{(A,B)\in M_3(\mathbb{F}_p)^2 | AB=BA\}|\leq \frac p{p-3}|G|\left(3p+\frac{p^3-p}3  +2p(p-1)+\frac{p(p^2-p)}2+\frac{p(p-1)(p-2)}6\right)
$$
$$
=\frac p{p-3} |G| (p^3+p^2+p).
$$
Combining these, we obtain that
$$
\frac{ |G|(p^3+p^2+p)}{(p^9)^2}\leq \frac{|\{(A,B)\in M_3(\mathbb{F}_p)^2 | AB=BA\}|}{(p^9)^2}\leq \frac p{p-3}\frac{|G|(p^3+p^3+p)}{(p^9)^2}.
$$
Since $|G|=(p^3-1)(p^3-p)(p^3-p^2)$, we have as $p\rightarrow\infty$, 
$$
\frac{|\{(A,B)\in M_3(\mathbb{F}_p)^2 | AB=BA\}|}{(p^9)^2}\sim \frac1{p^6}.
$$
