Find the cardinality of a subset of $GL_n( \mathbb F_p)$ 
Let $m,n \in \mathbb N$. Let $\mathbb F_p$ denote the prime field of characteristic $p$. Consider the set $$ X_m = \{A \in GL_n( {\mathbb F_p}): A^m=1 \}$$
  Compute the cardinality of $X_m$.

Its clear that $\vert X_m \vert < \infty$ since cardinality of $GL_n( \mathbb F_p)$ itself is $(p^n-1)(p^n-p)...(p^n-p^{n-1})$. Moreover, suppose $A \in X_m$ then $(x^m-1)$ kills $A$.
First I tried to understand the case when $m=p$. In this case if $A \in X_p$ then $(x^p-1)$ kills $A$ and since $x^p-1=(x-1)^p$ hence $(x-1)^n$ also kills $A$. Also, minimal polynomial of $A$ is of the form $(x-1)^k$ for some $k \leq p$. Any ideas to proceed further?   
 A: Suppose $m=qm'$ where $q$ is a power of $p$ and $m'$ is coprime to $p$. Then for $A\in X_m$ we can take the Jordan decomposition $A=us$, so  $u^ms^m=1$. By the uniqueness of the Jordan decomposition, $u^m=1$ and $s^m=1$. But $u$, being unipotent, has order that is a power of $p$, so $u^q=1$. Similarly, $s^{m'}=1$. Thus:
$$\# X_m=\sum_{s\in X_{m'}}\#\{u\in C_G(s):u^q=1\}.$$
Here, the eigenvalues of $s$ are multisets $Eig$ of Galois orbits of solutions to $X^m=1$. For each Galois orbit $\alpha\subset\mu_{m'}$, let $n_\alpha$ be the multiplicity of $\alpha$ in $Eig$, and let $q_\alpha:=p^{\deg(\alpha)}$ be the smallest power $q$ of $p$ such that $\alpha\subset\mathbb F_q$.
Now, we see that
$$\#\{u\in C_G(s):u^q=1\}=\prod_{\alpha\subset\mu_{m'}}\#\{u\in\mathrm{GL}_{n_\alpha}(\mathbb F_{q_\alpha}):u^q=1\}.$$
Moreover, given a semisimple element $s$, there are $\#(G/C_G(s))=\#\mathrm{GL}_n(\mathbb F_q)/\prod_\alpha\#\mathrm{GL}_{n_\alpha}(\mathbb F_{q^\alpha})$ semisimple elements conjugate to it, so in fact,
$$\# X_m=\sum_{Eig}\frac{\#\mathrm{GL}_n(\mathbb F_p)}{\prod_j\#\mathrm{GL}_{n_\alpha}(\mathbb F_{q_\alpha})}\prod_{\alpha\in Eig}\#\{u\in\mathrm{GL}_{n_\alpha}(\mathbb F_{q_\alpha}):u^q=1\}.$$
To better organize the formula, let
$$N_n^m(\mathbb F_q):=\#\{X\in\mathrm{GL}_n(\mathbb F_q):X^m=1\}/\#\mathrm{GL}_n(\mathbb F_q)$$ be the proportion of elements in $\mathrm{GL}_n(\mathbb F_q)$ with $X^m=1$. Now, the above formula can be re-written as
$$N_n^m(\mathbb F_p)=\sum_{Eig}\prod_{\alpha\in Eig} N_{n_\alpha}^q(\mathbb F_{q_\alpha}).$$
This can be re-phrased even further, by considering generating functions
$$f^m(\mathbb F_q)(T):=\sum_{n\ge0}N_n^m(\mathbb F_q)T^n.$$
Then, we have:
$$f^m(\mathbb F_p)(T)=\prod_{\alpha\subset\mu_{m'}}f^q(\mathbb F_{q_\alpha})(T^{n_\alpha}).$$
Now, the calculation of $\# X_m$ is essentially reduced to when $m$ is a power of a prime. When $m\ge p^n$ the set is exactly the set of nilpotent elements of $M_n(\mathbb F_p)$, which is $p^{n(n-1)}$, as mentioned in the comments. Otherwise, I do not know of any nice formulae for the size of such sets.

A particularly nice example is when $m|(p-1)$, in which we are looking for the coefficient of $T^n$ in the power series
$$\left(\sum_{k=0}^\infty\frac1{\#\mathrm{GL}_k(\mathbb F_p)}T^k\right)^m,$$
where $\#\mathrm{GL}_k(\mathbb F_p)=\prod_{i=0}^{k-1}(p^k-p^i)$ is some analog of the factorial. This function looks quite similar to the $q$-exponential, and is some hypergeometric series.
