why is $\sum\limits_{k=1}^{n} k^m$ a polynomial with degree $m+1$ in $n$ why is $\sum\limits_{k=1}^{n} k^m$ a polynomial with degree $m+1$  in $n$?
I know this is well-known. But how to prove it rigorously? Even mathematical induction does not seem so straight-forward.
Thanks.
 A: Qiaochu's answer is (as always?) right on the money.
I just wanted to mention that I also wrote up some notes about "discrete calculus" which takes the question of finding a closed form for the sums of $k$th powers as a jumping off point.  This came from a sophomore level undergraduate course I taught on introduction to proofs, in which we were doing suspiciously many induction proofs involving such closed forms, and I came to understand belatedly that the students found these confusing.  Not because they didn't understand the mechanism of induction (at least, not for the most part), but because they couldn't help but fixate on the nontrivial question: how do you know what closed form expression to put on the right hand side?  In the context of their exercises, the answer is: "You know because you are being given this information as part of the problem."  But of course that's a stupid answer: really, how do you know?
So I wrote up these notes on the subject of discrete calculus.  Or at least I started to: they are not quite complete and are a little more complicated than what I was expecting of the students, so I ended up not giving them out.  See especially Section 2, which explores the issue raised in the previous paragraph and also gives the extra little push (an innocuous recursive formula for the power sums) you need in order to prove by induction that $1^k + \ldots + n^k$ is a degree $k+1$ polynomial in $n$.  Note also that Qiaochu's nice linear algebra argument comes later, in Section 4 ("Linear Algebra of the Discrete Derivative") and especially Theorem 11.
As usual, any suggestions or corrections on these notes are most welcome, although (as usual) I certainly do not claim that there are few enough typos for it to be worth your time to point them out to me one by one.
A: Check Faulhaber's formula
A: Let $V$ be the space of all polynomials $f : \mathbb{N}_{\ge 0} \to F$ (where $F$ is a field of characteristic zero).  Define the forward difference operator $\Delta f(n) = f(n+1) - f(n)$.  It is not hard to see that the forward difference of a polynomial of degree $d$ is a polynomial of degree $d-1$, hence defines a linear operator $V_d \to V_{d-1}$ where $V_d$ is the space of polynomials of degree at most $d$.  Note that $\dim V_d = d+1$.  
We want to think of $\Delta$ as a discrete analogue of the derivative, so it is natural to define the corresponding discrete analogue of the integral $(\int f)(n) = \sum_{k=0}^{n-1} f(k)$.  But of course we need to prove that this actually sends polynomials to polynomials.  Since $(\int \Delta f)(n) = f(n) - f(0)$ (the "fundamental theorem of discrete calculus"), it suffices to show that the forward difference is surjective as a linear operator $V_d \to V_{d-1}$.
But by the "fundamental theorem," the image of the integral is precisely the subspace of $V_d$ of polynomials such that $f(0) = 0$, so the forward difference and integral define an isomorphism between $V_{d-1}$ and this subspace.  
More explicitly, you can observe that $\Delta$ is upper triangular in the standard basis, work by induction, or use the Newton basis $1, n, {n \choose 2}, {n \choose 3}, ...$ for the space of polynomials.  In this basis we have $\Delta {n \choose k} = {n \choose k-1}$, and now the result is really obvious.
The method of finite differences provides a fairly clean way to derive a formula for $\sum n^m$ for fixed $m$.  In fact, for any polynomial $f(n)$ we have the "discrete Taylor formula"
$$f(n) = \sum_{k \ge 0} \Delta^k f(0) {n \choose k}$$
and it's easy to compute the numbers $\Delta^k f(0)$ using a finite difference table and then to replace ${n \choose k}$ by ${n \choose k+1}$.  I wrote a blog post that explains this, but it's getting harder to find; I also explained it in my notes on generating functions.
A: You can set up a recursive formula for $\sum_{k=0}^n k^m $ by noting that
$$(n+1)^{m+1} = \sum_{k=0}^n (k+1)^{m+1}- \sum_{k=0}^n k^{m+1}$$
$$ = { m+1 \choose 1} \sum_{k=0}^n k^m
+ { m+1 \choose 2} \sum_{k=0}^n k^{m-1} + \cdots $$
by expanding the first summation on the RHS by the binomial theorem. Then shift all the other summations except $\sum_{k=0}^n k^m $ to the LHS.
So we get
$${ m+1 \choose 1} \sum_{k=0}^n k^m = (n+1)^{m+1} - { m+1 \choose 2} \sum_{k=0}^n k^{m-1}
- { m+1 \choose 3} \sum_{k=0}^n k^{m-2} + \cdots $$
A: The formula just drops right out if we use the Euler Maclaurin Summation Formula.
For $\displaystyle f(x) = x^m$ we have
$$ \sum_{k=0}^{n} f(k) = \int_{0}^{n} f(x)\ \text{d}x + \frac{n^m}{2} + \sum_{j=0}^{\infty} \frac{B_{2j}}{(2j)!} (f^{(2j-1)}(n) - f^{(2j-1)}(0))$$
Where $\displaystyle B_j$ are the Bernoulli numbers and $\displaystyle f^{(j)}(x)$ is the $\displaystyle j^{th}$ derivative of $\displaystyle f$.
Since $\displaystyle f(x)$ is polynomial, the terms in
$$ \sum_{j=0}^{\infty} \frac{B_{2j}}{(2j)!} (f^{(2j-1)}(n) - f^{(2j-1)}(0))$$
all are zero after a point ($\displaystyle 2j-1 \gt m$) and thus we get the formula for
$\displaystyle \sum_{k=0}^{m} k^m$ as a polynomial in $\displaystyle n$, with degree $m+1$.
A: $$ \sum_{k=1}^n k^m = \sum_{i=0}^{m-1} T(m,i) \binom{n+i+1}{m+1}$$
where $T(m,i)$ is an Eulerian number. (OEIS id:A008292)
Since $\binom{n+i+1}{m+1}$ is a polynomial of degree $m+1$ in $n$, so is $ \sum_{k=1}^n k^m $.
A: No one here has used induction via summation by parts?
Define $G_m(n) = \sum_{k=1}^n k^m$. We'll show that $G_m(n)$ is a polynomial of degree $m+1$. We know that $G_1$ is a polynomial of degree $2$ with leading coefficient of $1/2$. This is the basis step.
Step 1. For the inductive step, since $G_1(k)-G_1(k-1)=\frac{k(k+1)}{2}-\frac{(k-1)k}{2}=k$, summation by parts states that we can write $G_{m+1}(n)$ as
$$\sum_{k=1}^n [G_{1}(k)-G_1(k-1)] k^m  = [G_1(n+1)(n+1)^m - G_1(1)(1)^m] - \sum_{k=1}^n G_1(k)[(k+1)^m-k^m]$$
or if you prefer
$$G_{m+1}(n) = \left[\frac{(n+1)(n+2)}{2}(n+1)^m - 1 \right] - \sum_{k=1}^n \frac{k(k+1)}{2}[(k+1)^m-k^m].$$
Step 2. The first two terms are the boundary terms, as in normal integration by parts. Since $G_m(k)$ are polynomials of degree $m+1$ and $G_1$ has a leading coefficient of $1/2$, there are coefficients $a_j$ such that we can write
$$\sum_{k=1}^n G_1(k)[(k+1)^m-k^m] = \frac{1}{2} G_{m+1}(n) + \sum_{j=1}^m a_j G_j(n).$$
Step 3. Substituting and solving for $G_{m+1}(n)$ gives
$$\frac{3}{2} G_{m+1}(n)= \left[\frac{(n+1)(n+2)}{2}(n+1)^m - 1\right] - \sum_{j=1}^m a_j G_j(n).$$
Noting the orders of $G_j$, it is clear that $G_{m+1}$ has order $m+2$.
A: Let's define
$$S_p(n) = \sum_{k=0}^n k^p$$
Since $k^p=k\cdot k^{p-1}$, we can write
$$\begin{align}
S_p(n) = &1^{p-1} +\\
         &2^{p-1} + 2^{p-1} +\\
         &3^{p-1} + 3^{p-1} + 3^{p-1} +\\
         &\vdots\\
         &n^{p-1} + n^{p-1} + ... + n^{p-1}
\end{align}$$
We can add up this sum column by column. For $p=2$ and $p=3$ this has a nice interpretation as counting the number of cubes in stacks of 3d objects by going one planar cross-section at a time. If we do this, we get
$$\begin{align}
S_p(n) =& S_{p-1}(n) + \\
        &(S_{p-1}(n) - 1^{p-1}) + \\
        &(S_{p-1}(n) - 1^{p-1} - 2^{p-1}) +\\
        &\vdots \\
        &(S_{p-1}(n) - S_{p-1}(n-1))\\
       =&\sum_{k=0}^{n-1}\left(S_{p-1}(n) - S_{p-1}(k)\right)
\end{align}$$
It may be worthwhile to stop here and use this "add column by column" idea to derive $S_2(n)$ from $S_1(n)$, before reading further.
Now we need to assume for the sake of induction that, for all $k$ less than $p$, $S_k(n)$ is a polynomial in $n$ of degree $k+1$, with rational coefficients. This is of course true for $p=1$ as can be shown by Gauss's argument, and if you followed my advice in the last paragraph you'll know it's true for $p=2$ as well. Let's write
$$S_p(x) = a_{p+1}x^{p+1} + a_px^p + ... a_0$$
We can continue the calculation started above:
$$\begin{align}
S_p(n) &= \sum_{k=0}^{n-1}(S_{p-1}(n) - S_{p-1}(k))\\
       &= n S_{p-1}(n) - \sum_{k=0}^{n-1} S_{p-1}(k) \\
       &= n S_{p-1}(n) - \sum_{k=0}^{n-1} \sum_{i=0}^p a_i k^i \\
       &= n S_{p-1}(n) - \sum_{i=0}^{p} a_i S_i(n-1)
\end{align}$$
In the last step, we've taken the double summation and grouped the terms by exponent. You can see what we're beginning to do here - we're trying to express $S_p(n)$ in terms of $S_k(n)$ for $k<p$ so that it will follow by induction that if all of the latter are polynomials in $n$, so is $S_p(n)$. Right now the problem is that what we have are $S_k(n-1)$, not $S_k(n)$, but this is just fine since we have $S_k(n-1) = S_k(n) - n^k$. Thus, we may continue:
$$\begin{align}
S_p(n)&= n S_{p-1}(n) - \sum_{i=0}^{p} a_i (S_i(n) - n^i) \\
&= n S_{p-1}(n) - \sum_{i=0}^{p} a_i S_i(n) + \sum_{i=0}^p a_i n^i \\
&= (n + 1) S_{p-1}(n) - \sum_{i=0}^{p} a_i S_i(n)
\end{align}$$
Notice that the right hand side of this equation actually contains $S_p(n)$, so we have to "solve for" $S_p(n)$, obtaining:
$$S_p(n) = (1 - a_p)^{-1}\left((n + 1) S_{p-1}(n) - \sum_{i=0}^{p-1} a_i S_i(n)\right)$$
Notice that this is only valid if $a_p\neq -1$, so we need to amend our induction hypothesis: the $S_k(n)$ are not only polynomial in $n$ and of degree $k+1$, but they have leading coefficients not equal to $-1$. The inductive step immediately follows from our last equation. The leading coefficient of the polynomial in $n$ that is the right hand side is $\frac {a_p} {1 - a_p}$, and the function $\frac x {1 - x}$ sends positive reals to positive reals, so we can never get a leading coefficient of $-1$.
