The sum of same powers of all matrices modulo p The following is a problem from our department algebra competition for
students:

Non-question.
   An experimental-math geek was trying to raise all matrices $17\times17$
  over the field with 17 elements to the power of 100, sum the returns,
  and observe the result, when his computer broke. Help him.

Probably, the most of us can calculate the sum without use of computer
(alternatively, Russian readers can find a solution
here).
The problem seems to become much
harder when we replace 100 by, say, 80. More generally, my question is the
following:

Question.
  What is the sum
  $$
> \sum_{A\in M_p(\mathbb F_p)}A^k,
> $$
  where $p$ is prime and $k$ is a multiple of $p-1$?

It is easy to show that the sum is a scalar matrix, but which one?
Note that the coincidence of the matrix size and the characteristic
makes trace useless.
 A: The sum is zero for all $k<p^2-1$. 
Assume that $k$ is a multiple of $p-1$ and $k<p^2-1$. Divide all matrices into classes of the form $\{A,A+1,A+2,\dots,A+(p-1)\}$ . Summimg the $k$th powers over such a class yields
$$
 \sum_{s=0}^k A^{k-s}\binom ks \sum_{j=0}^{p-1} j^s
$$
where $0^0=1$. If $s=0$ or $s$ is not a multiple of $p-1$, the inner sum is 0 mod $p$. If $s$ is a multiple of $p-1$ and $s<k$, then the binomial coefficient is 0 mod $p$. (This is what breaks for $k=p^2-1$.) The only remaining term is the one $s=k$ and it does not depend on $A$. So the sum over a class is the same for all classes. The number of classes is a multiple of $p$, hence the result.
A: These are just several thoughts. but it seems that they show in particular that the answr is $0$ for $k< p^p-1$.
$\def\FF{{\mathbb F}}$
1. Firstly, denoting $d=\mathop{\rm lcm}(p-1,p^2-1,\dots,p^p-1)$ we see that $A^{pd+p}=A^p$ for every matrix $A$ (the order of a semisimple component divides $d$, and we need $p$ in order for a nilpotent component to vanish). Thus we may assume that $k  < pd+p$.
2. Your sum is equal to
$$
  S=\sum_{a_{11}\in\FF_p} \dots\sum_{a_{pp}\in\FF_p}
    \left(\sum_{i=1}^p\sum_{j=1}^p a_{ij}E_{ij}\right)^k,
$$
Now expand the inner brackets; we obtain
$$
  S=\sum_J \left(\sum_{a_{11}\in\FF_p} \dots\sum_{a_{pp}\in\FF_p} a^J\right)M_J,
$$
where the outer summation is taken over some multiindices $J$ with $|J|=k$, and $M_J$ are some matrices. The summation in brackets over $a_{ij}$ gives zero unless the exponent of $a_{ij}$ is nonzero and divisible by $p-1$. Thus the sum in the brackets vanishes unless all the coordinates of $J$ are nonzero and divisible by $p-1$. So, if $k < p^2(p-1)$ then $S=0$, as is in the case $k=80$ and $p=17$.
3. For the remaining case, we need to calculate $M_J$. Assume that $J=(j_{11},\dots,j_{pp})$ (all $j$'s are multiples of $p-1$). Consider a digraph $G_J$ with $\FF_p$ as the set of vertices and $j_{k\ell}$ edges from $k$ to $\ell$. Now, if $M_J=[m_{k\ell}]$ then $m_{kk}$ is the number of Eulerian paths starting at $k$ (multiplies by the sum in brackets which is $-1$), and all other entries are zeroes.
Now let us show that the number of such cycles is divisible by $p$. We will assume that $k=1$, so that the cycles start and end at $1$. Split each cycle into subcycles starting at $1$, ending at $1$ and not passing through $1$ any more. Correspond to each cycle all the cycles obtained by permutatons of subcycles; thus the set of all cycles is partitioned into such equivalence classes, and the number of elements in each class is a corresponding multinomial coefficient $\binom{s}{s_1,\dots,s_t}$ where $s_1,\dots,s_t$ are the numbers of occurences of different subcycles. This coefficient is divisible by $p$ unless in the $p$-adic notation there are no transitions in the addition $s_1+\dots+s_t=k$.
Notice that there are at least $p$ distinct subcycles --- at least one starting from each of the edges $1\to 1$, $1\to 2$, \dots, $1\to p$. Moreover, we may partition all the subcycles into such classes --- the number in each will be divisible by $p-1$. THus, we have $p$ nonzero numbers divisble by $p-1$, and there shpuld be no transitions in addition of these $p$ numbers; this may happen only if $k\geq (p-1)+p(p-1)+\dots+p^{p-1}(p-1)=p^p-1$. So, for $k< p^p-1$ we definitely have $S=0$. 
EDIT. Below in the comments, several improvements of this bound are shown.
A: [corrections applied, per Ilya and Anton]
Consider the formal series
$$
f(x) = \sum_{k=0}^\infty (\sum_A A^k) x^k
$$
where $A$ runs through $p \times p$ matrices.  It is equal to $\sum_A (I - Ax)^{-1}$, a rational function with values in scalar matrices.  Thus, for some $d$, if the first $d$ coefficients vanish all coefficients must vanish.
To determine $d$, we need to bound the degree of the numerator and denominator of $f(x)$.  Note by Cramer's rule, each $(I - Ax)^{-1}$ is a matrix of degree $(p-1)$ polynomials divided by the degree $p$ determinant of $(I - Ax)$.
There are $p^p$ different denominators, i.e. polynomials of the form 
$$
\det(I - Ax) = c_p x^p + c_{p-1} x^{p-1} + \cdots + c_1 x + 1
$$
The $c_i$ are unconstrained, they appear at least once in case $A$ is the companion matrix with $-c_p,-c_{p-1},\ldots$ along the last column.  So we can write the sum as
$$
\sum_{c_p,\ldots,c_1} \frac{\text{something of degree $p-1$}}{c_p x^p + \cdots + c_1 x + 1} 
$$
There are $(p-1)$ linear denominators, $p (p-1)$ quadratic denominators, $p^2 (p-1)$ cubic denominators and so on, so the product of all the $c_p x^p + \cdots + 1$ has degree
$$
d = (p-1) + 2p(p-1) + 3p^2(p-1) + \cdots + p p^{p-1}(p-1)
$$
i.e. $d = p^{p+1} - \frac{p^p - 1}{p - 1}$.  So $f(x)$ is a polynomial of degree less than $d$ divided by a polynomial of degree $d$.
Ilya's answer shows that the first $p^{p+1} - p$ Taylor coefficients of $f(x)$ vanish, and when $p \geq 3$ this is larger than $d$ so all the remaining Taylor coefficients are vanish as well.  When $p = 2$ Will points out the sums are not always zero, in fact for $p=2$
$$
f(x) = \frac{x^5}{(x+1)^2 (x^2 + x + 1)} I
$$
A: My answer was nonsense, sorry.
