How effectively Mathematica calculates Sum[Sum[a_k k^m, {m, 0, n}], {k, 1, p}] as a function of n?

Question

Suppose that I have a polynomial of order $n$ $$ f_n(k)=\sum_{m=0}^n a_k k^m, $$ where $k$ is an integer and $a_k$ are arbitrary real numbers. Now I want to use Mathemtica to calculate $$ \sum_{k=1}^p f_n(k). $$ Mathematica does it pretty fast even for large $n$.

Q0. How does Mathematica calculate these sums? I'm asking mainly about how the engine processes it. I suppose it takes $k^m$ (one by one $m=0,\dots,n$) and uses some known (analytical) procedures for summing the partial sums over $k$.

Q1. What is the scaling as a function of $n$ for a fixed $p$ Mathematica achieves? Some experimenting shows that is must be polynomial in time.

Q2. Independently on Mathematica, what is the computational complexity of calculating such sums as a function of $n$ (for a fixed $p$)?

EDIT: The coefficients of $k^m$ in $f_n(k)$ are set to one but they are otherwise arbitrary real numbers so $f_n(k)$ can't be in general simplified.

EDIT: I added coefficients $a_k$ to the first sum to make the question more clear.

Answer 1

I have no idea how Mathematica evaluates such sums. However, note that $f_n(k)=\frac{k^{n+1}-1}{k-1}$, meaning that

Q1. What is the scaling as a function of n for a fixed p Mathematica achieves? Some experimenting shows that is must be polynomial in time.

has the answer $O(1)$, provided that you force symbolic evaluation of $f_n(k)$.

Meanwhile, you could also ask what the scaling is as a function of $p$. It's natural to guess polynomial time, which you can verify by the following:

Sum[k^m, {m, 0, n}]
f[p_, n_] = Sum[%, {k, 2, p}]
DiscretePlot[Log[10, #] &@First@Timing@f[Round[10^k], 12], {k, 3, 6, 0.1}]

The slope of the line is 1, indicating linear time, or $O(p)$ execution.