A formula for the power sums: $1^n+2^n+\dotsc +k^n=\,$?

Is there explicit formula for the expression $1^n + 2^n + \dotsc + k^n\,$?

I know that for $n=1$ the explicit formula becomes $S=k(k+1)/2$ and for $n=3$ the formula becomes $S^2$. But what about general $n$?

I know there is a way using the Taylor expansion of $f(x)=1/(1-x)=1+x+x^2+\dotsc\;$, by differentiating it and then multiplying by $x$ and then differentiating again. Repeating this $n$ times, we get

$$\frac{d}{dx}(x\frac{d}{dx}(\dots x\frac{d}{dx}f(x))\dots )=1+2^nx^n+3^nx^n\dots.$$

Now do the same process but with the function $g(x)=x^{k+1}f(x)$. Then subtract them and we get $1+2^nx^n+\dots k^nx^n$. Because we have the explicit formulas $f(x)$ and $g(x)$ we can find the explicit formula by this process for arbitrary $n$. A big problem is that as $n$ grows, it is going take a lot of time finding the explicit formula. My question is therefore: are there other ways?

There is a result which does not depend on the Bernouli numbers or the Eulerian numbers. I found it in my first year at university, so you can be sure it involves nothing too complicated.

The method is that of discrete calculus, first I define $ \Delta f(x) = f(x+1)-f(x) $ and $ \Sigma f(x) = \sum\limits_{k=0}^{x-1} f(k) $ as discrete equivalents of differentiation and integration respectively.

The result will follow from a discrete equivalent of integration by parts, which I first proove:

$ f(x)g(x) = \Sigma f(x+1)g(x+1) -\Sigma f(x)g(x) $

$ \hspace{30 pt} = \Sigma f(x+1)g(x+1) -\Sigma f(x)g(x) + \Sigma f(x)g(x+1) - \Sigma f(x)g(x+1) $

$ \hspace{30 pt} = \Sigma (f(x+1)-f(x))g(x+1) + \Sigma f(x)(g(x+1)-g(x)) $

$ \hspace{30 pt} = \Sigma (\Delta f(x))g(x+1) + \Sigma f(x) \Delta g(x) $

$ \therefore \space \Sigma f(x) \Delta g(x) = f(x)g(x) - \Sigma (\Delta f(x))g(x+1) $

The expression we want a formula for will be a polynomial one degree higher than the power n, so we choose f(x) and g(x) with consideration that $ \deg(f(x)g(x)) $ should be $ n+1 $.

For a dash of simplicity, we can get $ f(x) \Delta g(x) = x^{n} $ by choosing $ f(x) = x^{n-1} $ and $ g(x) = \frac{1}{2} x(x-1) $.

For the other terms in the equivalence we need:

$ f(x)g(x) = \frac{1}{2} (x-1) x^n $

$ (\Delta f(x))g(x+1) = ((x+1)^{n-1}-x^{n-1})(\frac{1}{2} x(x+1)) = \frac{1}{2} x(x+1) \sum\limits_{m = 0}^{n-2} \binom{n-1}{m} x^m $

$ \hspace{57 pt} = \frac{1}{2} x(x+1)((n-1) x^{n-2} + \sum\limits_{m = 0}^{n-3} \binom{n-1}{m} x^m) $

$ \hspace{57 pt} = \frac{1}{2}(n-1) x^n + \frac{1}{2}(n-1) x^{n-1} + \frac{1}{2} \sum\limits_{m = 0}^{n-3} \binom{n-1}{m} (x+1) x^{m+1} $

Finally, plugging these in and rearranging:

$ \sum\limits_{k=0}^{x-1} k^n = \frac{1}{2} (x-1) x^n - \sum\limits_{k=0}^{x-1} \left( \frac{1}{2}(n-1) k^n + \frac{1}{2}(n-1) k^{n-1} + \frac{1}{2} \sum\limits_{m = 0}^{n-3} \binom{n-1}{m} (k+1) k^{m+1} \right) $

$ 2\sum\limits_{k=0}^{x-1} k^n = (x-1) x^n - (n-1) \sum\limits_{k=0}^{x-1} k^n - (n-1) \sum\limits_{k=0}^{x-1} k^{n-1} - \sum\limits_{k=0}^{x-1} \sum\limits_{m = 0}^{n-3} \binom{n-1}{m} (k+1) k^{m+1} $

$ (n+1) \sum\limits_{k=0}^{x-1} k^n = (x-1) x^n + (1-n) \sum\limits_{k=0}^{x-1} k^{n-1} - \sum\limits_{k=0}^{x-1} \sum\limits_{m = 0}^{n-3} \binom{n-1}{m} (k+1) k^{m+1} $

$ \sum\limits_{k=0}^{x-1} k^n = \frac{x-1}{n+1} x^n + \frac{1-n}{n+1} \sum\limits_{k=0}^{x-1} k^{n-1} - \sum\limits_{k=0}^{x-1} \sum\limits_{m = 0}^{n-3} \binom{n-1}{m} \frac{k+1}{n+1} k^{m+1} $

$ \therefore \sum\limits_{k=1}^{x} k^n = x^n + \sum\limits_{k=0}^{x-1} k^n = x^n + \frac{x-1}{n+1} x^n + \frac{1-n}{n+1} \sum\limits_{k=0}^{x-1} k^{n-1} - \sum\limits_{k=0}^{x-1} \sum\limits_{m = 0}^{n-3} \binom{n-1}{m} \frac{k+1}{n+1} k^{m+1} $

$ \hspace{31 pt} = \frac{x+n}{n+1} x^n + \frac{1-n}{n+1} \sum\limits_{k=1}^{x-1} k^{n-1} - \sum\limits_{k=1}^{x-1} \sum\limits_{m = 0}^{n-3} \binom{n-1}{m} \frac{k+1}{n+1} k^{m+1} $

Note that the power of k in the sums on the RHS does not exceed $ n-1 $.

Example with $ n = 3 $; the formula simplifies down to:

$$ \sum\limits_{k=1}^{x} k^3 = \frac{1}{4} (x^4 + 3x^3 -3 \sum\limits_{k=1}^{x-1} k^2 - \sum\limits_{k=1}^{x-1} k) $$

Which further simplifies to the correct polynomial. Only the partial sums for $k^2$ and $k$ need be known, and this formula will provide that of any higher degree.

A method which is more seldom used is that involving the Eulerian numbers. I found this solution myself by completely elementary means and "pattern-detection" only- so I liked it very much and I've made a small treatize about this. Unfortunately it is only in German, and since it is over 12 years old I don't want to translate it just now. However, the path to go and the development of the formulae go step by step, so it should be understandable/tractable and self-explanatory enough to see, how this can be done.

In effect, one finds a solution of the form

$$ P(m,n) = 1+2^m+3^m+...+n^m = e_{m,1}\binom{n+1}{m+1}+ e_{m,2}\binom{n+2}{m+1}+ ... +e_{m,m} \binom{n+m}{m+1} $$ with the Eulerian triangle $\{e_{r,c}\}_{r,c=0..\infty} $ $$ \begin{array} {} 1 \\ .&1 \\ .&1&1 \\ .&1&4&1 \\ .&1&11&11&1\\ ... \end{array} $$

We know that $$\sum_{i=1}^ki^{n+1}-(i-1)^{n+1}=k^{n+1}.$$

But the expression $i^{n+1}-(i-1)^{n+1}$ can be reduced to $$i^{n+1}-\left[\binom{n+1}{0}i^{n+1}-\binom{n+1}{1}i^{n}+...+(-1)^r\binom{n+1}{r}i^{n+1-r}...+(-1)^{n+1}\binom{n+1}{n+1}i^{0}\right]$$

which simplifies to $$\binom{n+1}{1}i^{n}+...+(-1)^{r-1}\binom{n+1}{r}i^{n+1-r}+...+(-1)^{n}\binom{n+1}{n+1}i^{0}.$$

Thus $$\sum_{i=1}^k\left[\binom{n+1}{1}i^{n}+...+(-1)^{r-1}\binom{n+1}{r}i^{n+1-r}+...+(-1)^{n}\binom{n+1}{n+1}i^{0}\right]=k^{n+1}.$$

Therefore $$(n+1)\sum_{i=1}^ki^{n}+\sum_{i=1}^k\left[-\binom{n+1}{2}i^{n-1}...+(-1)^{r-1}\binom{n+1}{r}i^{n+1-r}...+(-1)^{n}\right]=k^{n+1}.$$

We can then see that:$$(n+1)\sum_{i=1}^ki^{n}=k^{n+1}+\sum_{i=1}^k\left[\binom{n+1}{2}i^{n-1}...+(-1)^{r}\binom{n+1}{r}i^{n+1-r}...+(-1)^{n+1}\right].$$

Finally for all 2 $\leq r \leq n+1$, $$\sum_{i=1}^ki^{n}=\frac{1}{n+1}\left(k^{n+1}+\sum_{i=1}^k\left[\binom{n+1}{2}i^{n-1}...+(-1)^{r}\binom{n+1}{r}i^{n+1-r}...+(-1)^{n+1}\right]\right).$$

Note that this can be solved iteratively. E.g. once you find the value for n=1, you can use this for n=2, and n=3 and so on...

Some of the other answers say this in a more complicated way, but I've seen so much confusion around this question that I'd like to make it as clear as possible.

There is no good formula for $\sum_1^k i^n$ for the same reason that there is no good formula for $\int_0^x t(t-1)\ldots(t-n+1)dt$, namely that we are using the wrong basis for the space of polynomials.

For integrals, $t^n$ is a simpler choice (even an eigenbasis!): $\int_0^x t^n dt = \frac{1}{n+1} t^{n+1}$. If we want to compute the integral above, we will probably expand out and then use this formula.

With summation, or "discrete calculus", it is the other way around. $\sum_1^k i^n$ is messy, but $\sum_1^k i^{\underline{n}} = \sum_1^k i(i-1)\ldots(i-n+1)$ is very simple, just $\frac{1}{n+1} (k+1)^{\underline{n+1}} = \frac{1}{n+1} (k+1)k(k-1)\ldots (k-n+1)$. (This doesn't give us an eigenbasis, but that's because we should really be summing up to $k-1$.)

I prefer to write this instead using binomial coefficients, ${i\choose n} = \frac{1}{n!}i^{\underline{n}}$. Then the formula is $\sum_1^k {i\choose n} = {{k+1}\choose{n+1}}$. This is easy to prove with a combinatorial argument or a number of other methods (see my post here).

But any formula for $\sum_1^k i^n$ will involve some complexity, and most likely be equivalent to expanding in terms of the "correct" basis. Even the nice-looking Eulerian numbers formula has two layers of complexity, one in the outer sum and one in the combinatorics of the numbers themselves.

At the bottom, this is very close to the other approaches, but the details are somewhat different, and it adds a nice perspective to the question:

Consider the function $$ B(t,x)=\frac{te^{tx}}{e^t-1}. $$ Expanding as a power series o=in $t$, we see that $$ B(t,x)=\sum_{m=0}^\infty B_m(x)\frac {t^m}{m!} $$ for some $B_m(x)$. (We can treat this formally, or check that it has a positive radius of convergence, and argue there.)

Note that $B(t,x+1)=te^{xt}+B(t,x)$, so $$ \sum_{m=0}^\infty(B_m(x+1)-B_m(x))\frac{t^m}{m!}=te^{xt}=\sum_{m=0}^\infty\frac{mx^{m-1}t^m}{m!}, $$ or $B_m(x+1)-B_m(x)=mx^{m-1}$.

From this it follows easily that each $B_m(x)$ is a polynomial in $x$ of degree $m$, with rational coefficients; these are the Bernoulli polynomials, and we have $$ \sum_{i=0}^n (B_{m+1}(i+1)-B_{m+1}(i))=\sum_{i=0}^n(m+1)i^m, $$ so $$ \sum_{i=0}^n i^m = \frac{B_{m+1}(n+1)-B_{m+1}(0)}{m+1}. $$

From a personal perspective, let me add that this version of the formula led to a nice result, discussed here: For any integer $k>2$ there are only finitely many primes of the form $p=kn+1$ such that $1^k,2^k,\dots,n^k$ are distinct modulo $p$. (This is false for $k=1,2$, and related to a nice open problem.)