Methods to compute $\sum_{k=1}^nk^p$ without Faulhaber's formula

Question

As far as every question I've seen concerning "what is $\sum_{k=1}^nk^p$" is always answered with "Faulhaber's formula" and that is just about the only answer. In an attempt to make more interesting answers, I ask that this question concern the problem of "Methods to compute $\sum_{k=1}^nk^p$ without Faulhaber's formula for fixed $p\in\mathbb N$". I've even checked this post of common questions without finding what I want.

Rule #1: Any method to compute the sum in question for arbitrary $p$ is good, either recursively or in some manner that is not in itself a closed form solution. Even algorithms will suffice.

Rule #2: I don't want answers confined to "only some values of $p$". (A good challenge I have on the side is a generalized geometric proof, as that I have not yet seen)

Exception: If your answer does not generalize to arbitrary $p$, but it still generalizes to an infinite amount of special $p$'s, that is acceptable.

Preferably, the method is to be easily applied, unique, and interesting.

To start us off, I have given my answer below and I hope you all enjoy.

Answer 1

Feel free to skip to the highlighted parts and the ending to see the formula in action.

---

Suppose we had a continuous and differentiable function that satisfied the following equation:

$$f(x,p)=f(x-1,p)+x^p,\quad f(0,p)=0$$

Differentiating with respect to $x$, we get

$$f'(x,p)=f'(x-1,p)+px^{p-1}$$

Now notice that when $x\in\mathbb N$

$$f(x,p)=\sum_{k=1}^xk^p$$

$$f'(x,p)=f'(0,p)+p\sum_{k=1}^xk^{p-1}=f'(0,p)+pf(x,p-1)$$

Integrating both sides from $0$ to $x$, we have

$$f(x,p)=\int_0^xf'(t,p)dt=\int_0^xf'(0,p)+pf(t,p-1)dt=xf'(0,p)+p\int_0^xf(t,p-1)dt$$

When $x=1$, it is easy enough to see that

$$a_p=f'(0,p)=1-p\int_0^1f(t,p-1)dt$$

Combining all of this, we have

$$f(x,p)=a_px+p\int_0^xf(t,p-1)dt$$

And with $p=0$, it is trivial to see that

$$a_0=1\implies f(x,0)=x$$

---

Further,

$$a_1=1-\int_0^1t\ dt=\frac12$$

$$f(x,1)=\frac12x+\int_0^xt\ dt=\frac12x+\frac12x^2$$

---

$$a_2=1-2\int_0^1\frac12t+\frac12t^2dt=\frac16$$

$$f(x,2)=\frac16x+\frac12x^2+\frac13x^3$$

Indeed, these are the solutions to the sum in question found by a recursive formula involving integrals.

---

This method is described here.

Answer 2

With $[z^n]$ denoting the coefficient of $z^n$ in a series and $D_z:=\frac{d}{dz}$ we obtain \begin{align*} \sum_{k=1}^nk^p=[z^n]\frac{1}{1-z}(zD_z)^p\frac{1}{1-z}\qquad\qquad\qquad p\in\mathbb{N} \end{align*}

See method 1 in this answer which derives this formula based upon generating functions together with a small example ($n=2$).

Another variation:

The sum of the $p$-th powers of numbers $1$ to $n$ is given by \begin{align*} \sum_{k=1}^nk^p=\sum_{k=1}^p{p\brace k}\frac{(n+1)^{\underline{k+1}}}{k+1}\qquad\qquad\qquad p\in\mathbb{N} \end{align*}

See method 2 in this answer together with a small example ($n=2$).

Here we use the Stirling Numbers of the second kind ${n\brace k}$ and Don Knuths falling factorial power notation: $n^{\underline{k}}=\frac{n!}{(n-k)!}$.

Answer 3

By the binomial theorem, $$(x+1)^{n+1}=\sum_{h=0}^{n+1} {n+1 \choose h}x^h$$ $$(x+1)^{n+1}-x^{n+1}=\sum_{h=0}^n {n+1 \choose h}x^h$$ Sum this equality for $x=0,1\dotsb k$ $$\sum_{x=1}^k((x+1)^{n+1}-x^{n+1})=(k+1)^{n+1}-1=\sum_{x=1}^k\sum_{h=0}^n {n+1 \choose h}x^h=\sum_{h=0}^n{n+1 \choose h}\sum_{x=1}^kx^h=(n+1)\sum_{x=1}^kx^h+\sum_{h=0}^{n-1}{n+1 \choose h}\sum_{x=1}^kx^h$$

Which means $$(n+1)\sum_{x=1}^kx^n=(k+1)^{n+1}-1-\sum_{h=0}^{n-1}{n+1 \choose h}\sum_{x=1}^kx^h$$

So you can find the sum of the $n$th powers if you have all the previous ones. The base case is $$\sum_{x=1}^kx^0=k$$

Then $$\sum_{x=1}^kx^1=\frac{1}{2}\left((k+1)^2-1-{2 \choose 0}k\right)=\frac{k^2+k}{2}$$

$$\sum_{x=1}^kx^2=\frac{1}{3}\left((k+1)^3-1-{3 \choose 0} k - {3 \choose 1} \frac{k^2+k}{2}\right)=\frac{k^3}{3}+\frac{k^2}{2}+\frac{k}{6}$$ Etcetera.

Answer 4

$(1)$

You are using the method of differences of Discrete Mathematics. For sums exist an equivalent to the integration. One of the pioneers of this method was Newton with the so-called Newtonian interpolation but in our days it’s written more elegant. You can read about this method in books of Concrete Mathematics, e.g. :

Ronald L. Graham, Donald E. Knuth, Oren Patashnik $ \enspace / \enspace $ Concrete Mathematics, Second Edition (2009) $ \enspace / \enspace $ Addison-Wesley Publishing Company, 1994 $ \enspace / \enspace $ http://www-cs-faculty.stanford.edu/~uno/gkp.html

Now a bit more general case for calculating the sum $\enspace \sum\limits_{k=1}^n k^p$ :

Be $f(x)$ any polynomial, $I$ the Identity operator with $If(x):=f(x)$,

$E$ the Displacement operator with $Ef(x):=f(x+1)$ and

$\Delta$ the Difference operator with $\Delta:=E-I$ .

We get $\enspace \Delta f(x)=(E-I)f(x)=Ef(x)-If(x)=f(x+1)-f(x)\enspace $

and with $\enspace \Delta^{n+1}:=\Delta(\Delta^n)\enspace $ the non-trivial formula $$\Delta^n f(x)=(E-I)^n f(x)= \sum\limits_{k=0}^n (-1)^{n-k}\binom{n}{k}E^k f(x) =\sum\limits_{k=0}^n (-1)^{n-k}\binom{n}{k}f(x+k)\enspace .$$

Be $\enspace\displaystyle m,n\in\mathbb{N}_0\enspace $ and $\enspace\displaystyle f_m(x):= \sum\limits_{k=0}^m b_k x^\underline{k}\in\mathbb{R}[x]\enspace $ a polynomial of degree $\enspace m$

with $\enspace\displaystyle x^\underline{n} :=\prod\limits_{v=1}^n (x-v+1)\enspace $ and therefore $\enspace\displaystyle \Delta^k x^\underline{n} =n^\underline{k} x^\underline{n-k} $ .

Then it follows $$\sum\limits_{j=0}^n f_m(j)= \sum\limits_{k=0}^m \binom{n+1}{k+1} \sum\limits_{v=0}^k (-1)^{k-v}\binom{k}{v}f_m(v) \enspace .$$

Proof:

With $\enspace\displaystyle \Delta^v f_m(x)=\sum\limits_{k=v}^m b_v k^\underline{v} x^\underline{k-v} \enspace $ and $\enspace x:=0\enspace $ follows $\enspace\displaystyle b_k=\frac{1}{k!} \Delta^k f_m(0) \enspace $ and therefore $\enspace\displaystyle f_m(x)= \sum\limits_{k=0}^m \frac{ x^\underline{k} }{k!} \Delta^k f_m(0) = \sum\limits_{k=0}^m \binom{x}{k}\sum\limits_{v=0}^k (-1)^{k-v}\binom{k}{v}f_m(v) $ .

Summation from $x=0$ to $n$ leads to $\enspace\displaystyle \sum\limits_{x=0}^n \binom{x}{k}= \binom{n+1}{k+1}\enspace $ and therefore to the confirmation of the assertion above.

With $\enspace m:=p\enspace $ and $\enspace f_p(x):=x^p\enspace $ you get a well-known formula for your sum: $$\sum\limits_{j=0}^n j^p= \sum\limits_{k=0}^p \binom{n+1}{k+1}{p\brace k}k!$$

(You know: ${p\brace k}$ is called the Stirling number of the second kind and gives you a better overview to compare with the other answers)

---

$(2)$

A complete other technique is the Euler-Maclaurin formula for summation approximation. It gives your sum a representation with the Bernoulli numbers.

But in the case of the concerned sum it's much better to calculate directly with the Bernoulli polynomials $B_k(x)$ defined by $\enspace\displaystyle \sum\limits_{k=0}^\infty \frac{B_k(x)}{k!}z^k:=\frac{ze^{xz}}{e^z-1}\enspace$ (generating function).

You have done this in your post, please look at my comment there.

It's $\enspace\displaystyle B_k(x)=\sum\limits_{v=0}^k \binom{k}{v}B_v x^{k-v}\enspace$ with $\enspace\displaystyle B_k=-\frac{1}{k+1}\sum\limits_{v=0}^{k-1}\binom{k+1}{v}B_v\enspace$, $\enspace k\in\mathbb{N}\enspace $ and $\enspace B_0=1\,$ .

It follows $$\sum\limits_{k=0}^\infty \frac{B_k(x+1)-B_k(x)}{k!}z^k=\frac{z}{e^z-1}(e^{(x+1)z}-e^{xz})=ze^{xz}=\sum\limits_{k=0}^\infty \frac{x^kz^{k+1}}{k!}$$ and therefore $\enspace B_k(x+1)-B_k(x)=kx^{k-1}$.

With $\enspace k\to p+1\enspace $ and $\enspace x\to k\enspace $ we get $$\sum\limits_{k=1}^n k^p =\sum\limits_{k=1}^n \frac{B_{p+1}(k+1)-B_{p+1}(k)}{p+1}= \frac{B_{p+1}(n+1)-B_{p+1}(1)}{p+1}$$ which is $\enspace \int\limits_1^{n+1}B_p(x)dx \enspace$ too.

---

$(3)\enspace $ More general considerations can be found in The sum of fractional powers $\sum\limits_{k=1}^x k^t$. .

Answer 5

The following gives a mostly impractical means of computing this sum.

Consider the Discrete Uniform Distribution with support on $\left\{1,2,\ldots,n\right\}$. Let $X$ be a random variable with this distribution.

Then $$E(X^p) = \frac{\sum_{k=1}^{n}k^p}{n}$$

The Moment Generating Function for this distribution is

$$E(e^{tX}) = \frac{e^{t} - e^{(n+1)t}}{n(1-e^t)}$$

then

$$\sum_{k=1}^{n}k^p = nE(X^p) = \frac{\operatorname{d}^p}{\operatorname{d}t^p}\left|_{t=0}\frac{e^{t} - e^{(n+1)t}}{1-e^t}\right.$$

If you would like to generalize this method to nonnegative real $p$, you may use Fractional Differentiation in the above formula, though this is even more impractical.

I'm pretty sure Faulhaber's formula can be derived with this method by using the General Leibniz Rule and Faà di Bruno's formula.

Let $f(t) = e^{t} - e^{(n+1)t}$, $g(t) = 1 - e^t$, and $u(t) = \frac{1}{g(t)}$.

$$\frac{\operatorname{d}^p}{\operatorname{dt^p}}\left(\frac{f(t)}{g(t)}\right) = \frac{\operatorname{d}^p}{\operatorname{dt^p}}\left(f(t)u(t)\right) = \sum_{k=0}^p{p \choose k}f^{(k)}(t)u^{(p-k)}(t)$$

Now $$f^{(k)}(t) = n^ke^{nt} - 2^{k}e^{2t}$$

While by Faà di Bruno's formula

$$u^{(p-k)}(t) = \sum_{r=0}^{p-k}(-1)^r\frac{r!}{g(t)^{r+1}}B_{p-k,r}\left(g^{(1)}(t),g^{(2)}(t),\ldots,g^{(p-r+1)}(t)\right)$$

where the $B_{(n-p,r)}$ are Ordinary Bell Polynomials.

Now to get our answer, substitute in $t = 0$.

Note that

$$B_{p-k,r}\left(g^{(1)}(0),\ldots,g^{(p-k-r+1)}(0)\right) = B_{p-k,r}\left(1,1,\ldots,1\right) = \left\{\begin{array}{c}p-k\\r\end{array}\right\}$$

Where $\left\{\begin{array}{c}p-k\\r\end{array}\right\}$ is a Stirling Number of the Second Kind.

Giving us the formula

$$\sum_{k=1}^{n}k^p = \lim_{t\rightarrow 0} \sum_{k=0}^{p}\sum_{r=0}^{p-k}(-1)^r r! {p \choose k}\left\{\begin{array}{c}p-k\\r\end{array}\right\}\left( \frac{e^{t} - (n+1)^ke^{(n+1)t}}{(1 - e^t)^{r+1}}\right)$$

Let's make the substitution $z = 1 - e^t)$ then we now seek

$$\lim_{z\rightarrow 0} \sum_{k=0}^{p}\sum_{r=0}^{p-k}(-1)^r r! {p \choose k}\left\{\begin{array}{c}p-k\\r\end{array}\right\} \frac{1-z - (n+1)^k(1-z)^{n+1}}{z^{r+1}}$$

Interchanging the order of summation and finding a common denominator gives us

$$\lim_{z\rightarrow 0} \frac{1}{z^{p+1}}\sum_{r=0}^{p}\sum_{k=0}^{p-r}(-1)^r r! {p \choose k}\left\{\begin{array}{c}p-k\\r\end{array}\right\} \left[1-z - (n+1)^k(1-z)^{n+1}\right]z^{p-r}$$

Since we know the limit must exist, since

$$\frac{e^t - e^{(n+1)t}}{1-e^t} = e^t + e^{2t} + \cdots + e^{nt}$$

which can be differentiated infinitely many times, we can simply find the coefficient of $z^{p+1}$ in the numerator to compute the limit.

By using the binomial expansion of $(1-z)^{n+1}$ and splitting up the sum, the numerator can be expressed as

$$\sum_{r=0}^{p}\sum_{k=0}^{p-r}(-1)^{r}r!{p \choose k}\left\{\begin{array}{c}p-k\\r\end{array}\right\}z^{p-r} + \sum_{r=0}^{p}\sum_{k=0}^{p-r}(-1)^{r}r!{p \choose k}\left\{\begin{array}{c}p-k\\r\end{array}\right\}z^{p-r+1} \:\:- $$ $$\sum_{r=0}^{p}\sum_{k=0}^{p-r}\sum_{j=0}^{n+1}(-1)^{r}r!{p \choose k}\left\{\begin{array}{c}p-k\\r\end{array}\right\}z^{p+j-r}$$

The first sum contributes nothing to the coefficient of $z^{p+1}$ since $p-r$ never equals $p+1$ when $0 \leq r \leq p$. The second sum can only contribute when $r = 0$, but in this case $\left\{\begin{array}{c}p-k\\0\end{array}\right\} = 0$ and we see that this sum does not contribute.

The third sum contributes when $j = r+1$, after some manipulations yielding

$$\sum_{r=0}^{p}\sum_{k=0}^{p-r}r!{p \choose k}\left\{\begin{array}{c}p-k\\r\end{array}\right\}{n+1 \choose r+1}$$

Now make the substitution $i = k+r$. The sum can then be brought to the form

$$\sum_{r=0}^{p}\sum_{i=0}^{p-r}r!{p \choose p+r-i}\left\{\begin{array}{c}p+r-i\\r\end{array}\right\}{n+1 \choose r+1}$$

Using the identity $$\left\{\begin{array}{c}n+1\\k+1\end{array}\right\} = \sum_{j=k}^{n}{n \choose j}\left\{\begin{array}{c}j\\k\end{array}\right\}$$

This can be brought to the form

$$\sum_{r=0}^{p}r!{n+1 \choose r+1}\left\{\begin{array}{c}p+1\\r+1\end{array}\right\}$$

which very nearly agrees with a formula given in this answer. I seem to have made an off by one error of some kind. I will try to fix it later.

I'll try to produce an example using fractional derivatives tomorrow.

Answer 6

We'll use Lagrange interpolation theorem.

We generalize the formula by letting $S_p:{\mathbb Z}\to {\mathbb Z}$ be a function satisfying:

1. $S_p(n+1) -S_p(n) = (n+1)^p,~n\in {\mathbb Z}$.
2. $S_p(0) =0$.

It follows that $S_p$ is a polynomial of degree $p+1$ (I'm skipping the proof to get to the main course).

By Lagrange interpolation theorem we can express this polynomial explicitly through the values it assigns to the $p+2$ points $\{-1,0,1,\dots,p\}$. Note that by 1. and 2. $S_p(0)=S_p(-1)=0$. The formula from the interpolation theroem is:

$$S_p(n) = \sum_{j=1}^p S_p(j)\times\prod_{k \in \{-1,0,\dots,p\}-\{j\}}\frac{n-k}{j-k},$$

which can be also expressed as

$$S_p(n) = \left( \prod_{k=-1}^p (n-k)\right) \sum_{j=1}^p \left (\prod_{k\in \{-1,\dots,p\}-\{j\}} \frac{1}{k-j} \times \frac{S_p(j)}{n-j} \right).$$

That's it.

Let's do some examples. For $S_1$ and $S_2$ we use the first variant, while for $S_3$ we use the second variant, skipping the intermediate calculations.

$$S_1(n) = 1 \frac{n+1}{1+1}\times \frac{n-0}{1-0}=\frac{n(n+1)}{2}.$$

\begin{align*} S_2(n) &= 1\times \frac{n+1}{2}\times \frac{n}{1}\times \frac{n-2}{1-2} \\ & + (1+4) \times \frac{n+1}{3}\times \frac{n}{2}\times \frac{n-1}{1}\\ & = \frac{(n+1)n}{6} \left (-3 (n-2) + 5 (n-1) \right)\\ & = \frac{(n+1)n (2n+1)}{6}. \end{align*}

\begin{align*} S_3(n) = (n+1)n (n-1)(n-2)(n-3) \times \left( \frac{1}{4(n-1)} -\frac{3}{2(n-2)}+\frac{3}{2(n-3)} \right), \end{align*} noting that to get to the nice formula for $S_3(n)$ as $S_1(n)^2$, we need to do more algebra.

Answer 7

Very elementary method: knowing that $S(n) = \sum_{k=1}^n k^p = a_{p+1}n^{p+1} + \cdots + a_1 n + a_0$, you can calculate the coefficients using limits: $$a_{p+1} = \lim_{n\to\infty}\frac{S(n)}{n^{p+1}} = \lim_{n\to\infty}\frac{1^p+\cdots+n^p}{n^{p+1}} = \lim_{n\to\infty}\frac{(n+1)^p}{(n+1)^{p+1}-n^{p+1}} = \cdots = \frac1{p+1}.$$ (Cesàro-Stolz used in the third =)

You can continue with $$a_p = \lim_{n\to\infty}\frac{S(n)-a_{p+1}n^{p+1}}{n^p} = \cdots$$ $$\cdots$$

Answer 8

The three methods I usually use to compute the formulas for $$ \sum_{k=1}^nk^m $$ use the Binomial Theorem, the Euler-Maclaurin Sum Formula, and the Hockey-Stick Identity.

---

Binomial Theorem

This recursive method requires one to compute the formulas for the sum of all $j^\text{th}$ powers for $j\lt m$ before computing the formula for the the sum of the $m^\text{th}$ powers.

Summing the binomial identity $$ k^{m+1}-(k-1)^{m+1}=\sum_{j=0}^m\binom{m+1}{j}(-1)^{m-j}k^j $$ gives $$ \begin{align} n^{m+1} &=\sum_{j=0}^m\binom{m+1}{j}(-1)^{m-j}\sum_{k=1}^nk^j\\ &=(m+1)\sum_{k=1}^nk^m+(-1)^mn+(m+1)\sum_{j=1}^{m-1}\binom{m}{j-1}\frac{(-1)^{m-j}}j\sum_{k=1}^nk^j \end{align} $$ Therefore, $$ \bbox[5px,border:2px solid #C0A000]{\sum_{k=1}^nk^m=\frac{n^{m+1}+(-1)^{m+1}n}{m+1}-\sum_{j=1}^{m-1}\binom{m}{j-1}\frac{(-1)^{m-j}}j\sum_{k=1}^nk^j} $$

---

Euler-Maclaurin Sum Formula

The Euler-Maclaurin Sum Formula is $$ \sum_{k=1}^nf(k)=\int f(n)\,\mathrm{d}n+\frac12\,f(n)+\sum_{j=1}^\infty\frac{B_{2j}}{(2j)!}\,f^{(2j-1)}(n) $$ Setting $f(k)=k^m$, we get $$ \bbox[5px,border:2px solid #C0A000]{\sum_{k=1}^nk^m=\frac{n^{m+1}}{m+1}+\frac{n^m}2+\frac1{m+1}\sum_{j=1}^\infty B_{2j}\binom{m+1}{2j}\,n^{m-2j+1}} $$ which turns out to be Faulhaber's Formula.

---

Hockey-Stick Identity

We can write the the monomial $k^m$ as a linear combination of the binomial coefficients $\binom{k}{j}$, which can be considered as a polynomial in $k$ of degree $j$: $$ \newcommand{\stirtwo}[2]{\left\{#1\atop#2\right\}} k^m=\sum_{j=0}^m\binom{k}{j}\stirtwo{m}{j}j! $$ where $\stirtwo{m}{j}$ are Stirling Numbers of the Second Kind.

Using the Hockey-Stick Identity, we get $$ \bbox[5px,border:2px solid #C0A000]{\sum_{k=1}^nk^m=\sum_{j=0}^m\binom{n+1}{j+1}\stirtwo{m}{j}j!} $$

Answer 9

A purely combinatorial approach:

Notice that $\sum\limits_{k=1}^n k^p$ is equal to the number of $p$-tuples $(a_1,a_2,\dots,a_p)$ such that $a_1\geq a_j$ for all $j$ and every term is among $\{1,2,\dots,n\}$.

We count them as follows:

How many such $p$-tuples contain exactly $j$ different terms? There are $p \brace j$ ways to split the $p$ terms into $j$ unlabelled groups and $\binom{n}{j}$ ways to select the $j$ different numbers, after this there are $(j-1)!$ ways to assign the $j$ numbers among the groups (since the group containing $a_1$ must be assigned the largest number).

We conclude $\sum\limits_{k=1}^nk^p=\sum\limits_{j=1}^p {p\brace j}\binom{n}{j} (p-1)!$.

Answer 10

We all know the geometric series, and from it, we can derive the following:

$$\sum_{k=1}^ne^{kx}=\frac{1-e^{(n+1)x}}{1-e^x}-1$$

Upon differentiating $p$ times and taking the limit as $x\to0$, we get

$$\sum_{k=1}^nk^p=\lim_{x\to0}\frac{d^p}{dx^p}\frac{1-e^{(n+1)x}}{1-e^x}$$

If one uses the Riemann–Liouville integral, then an integral form for negative $p$ can be created.

$$\sum_{k=1}^n\frac1{k^p}=\frac1{\Gamma(p)}\int_0^\infty\frac{1-e^{-nx}}{e^x-1}x^{p-1}~{\rm d}x$$

Answer 11

A "method",

For a closed formula and proof of the method used in this answer Formula for $1^k+2^k+3^k...n^k$ for $n,k \in \mathbb{N}$.

With the ideas used in answer it is simple to generate a formula for the sum of any power using this method,

The method given if power is $n$,

1. List $a_0=0$ and the first $n$ terms
2. Compute Finite differences
3. Construct polynomial in terms of binomial coefficients
4. Apply Hockey Stick Identify

This is best explained with an example, take $n=3$, then we do step one,

$$0,1^3,2^3,3^3$$

Now step $2$, take differences of term and previous term.

$$1^3-0,2^3-1^3,3^3-2^3$$

$$(2^3-1^3)-(1^3-0),(3^3-2^3)-(2^3-1^3)$$

$$(3^3-2^3)-(2^3-1^3)-((2^3-1^3)-(1^3-0))$$

Note the first terms of each sequence including the original is $0,1,6,6$ respectively this will correspond with step $3$.

Step $3$,

$$x^3=0{x \choose 0}+1{x \choose 1}+6{x \choose 2}+6{x \choose 3}$$

Step $4$,

$$\sum_{x=1}^{n} x^3=1{n+1 \choose 1+1}+6{n+1 \choose 2+1}+6{n+1 \choose 3+1}$$

Or if you prefer,

$$\sum_{x=1}^{n} x^3=0\frac{n+1}{1!}+1\frac{(n+1)(n)}{2!}+6\frac{(n+1)(n)(n-1)}{3!}+6\frac{(n+1)(n)(n-1)(n-2)}{4!}$$

Answer 12

This is a simple presentation of the Newton series that has been mentioned in user90369's answer.

---

The forward difference operator $Δ$ defined as: $ \def\nn{\mathbb{N}} \def\zz{\mathbb{Z}} \def\lfrac#1#2{{\large\frac{#1}{#2}}} \def\lbinom#1#2{{\large\binom{#1}{#2}}} $

$Δ = ( \text{function $f$ on $\zz$} \mapsto ( \text{int $n$} \mapsto f(n+1) - f(n) ) )$

Namely, for any function $f$ on $\zz$ and $n \in \zz$, we have $Δ(f)(n) = f(n+1) - f(n)$.

If you think of the functions as sequences (infinite in both directions), then taking the forward difference means replacing each term with the value of the next term minus itself. For example if you repeatedly take the forward difference of the sequence of cubes:

...,-27,-8,-1, 0, 1, 8,27,...
..., 19, 7, 1, 1, 7,19,37,...
...,-12,-6, 0, 6,12,18,24,...
...,  6, 6, 6, 6, 6, 6, 6,...
...,  0, 0, 0, 0, 0, 0, 0,...
...,  0, 0, 0, 0, 0, 0, 0,...

This powerful abstraction makes it easy to get a lot of things. For instance, the numbers obtained here can be easily used to obtain the general formula for sum of cubes, as you desire.

General method for indefinite summation

The key is that:

$Δ\left( \text{int $n$} \mapsto \lbinom{n}{k+1} \right) = \left( \text{int $n$} \mapsto \lbinom{n}{k} \right)$ for any $k \in \zz$.

This is to be expected because it follows directly from extending Pascal's triangle naturally, namely if we define $\lbinom{n}{k}$ by the recurrence:

$\lbinom{n}{0} = 1$ for any $n \in \zz$.

$\lbinom{0,k+1}{0} = 0$ for any $k \in \nn$.

$\lbinom{n+1}{k+1} = \lbinom{n}{k+1} + \lbinom{n}{k}$ for any $k \in \nn$ and $n \in \zz$.

Now consider any function $f$ on $\zz$ such that $f(n) = \sum_{k=0}^d a_k \lbinom{n}{k}$ for any $n \in \zz$. Then we have for any $m \in \nn_{\le d}$:

$Δ^m(f)(n) = \sum_{k=0}^{d-m} a_{k+m} \lbinom{n}{k}$ for any $n \in \zz$.

And hence:

$Δ^m(f)(0) = a_m$.

Which immediately gives Newton's series:

$f(n) = \sum_{k=0}^d Δ^k(f)(0) \lbinom{n}{k}$ for any $n \in \zz$.

From a high-level perspective, this is the discrete version of the Taylor series.

This works for any polynomial function $f$ on $\zz$, since $D^k(f)$ is the zero function once $k$ is larger than the degree of $f$, so we can use it to immediately find the series for $(\text{int n} \mapsto n^3)$, and then just take the anti-difference by shifting the coefficients of the series the other way. The undetermined constant that appears will drop out once we perform a definite sum like if we want the sum of the first $m$ cubes.

---

Sum of $p$ powers

For example if we want $\sum_{k=1}^{n-1} k^3$ we first find the iterated forward differences of the sequence of cubes $( n^3 )_{n \in \zz}$:

..., 0, 1, 8,27,64,...
..., 1, 7,19,37,...
..., 6,12,18,...
..., 6, 6,...
..., 0,...

So we immediately get $n^3 = 0 \binom{n}{0} + 1 \binom{n}{1} + 6 \binom{n}{2} + 6 \binom{n}{3}$ and hence $\sum_{k=0}^{n-1} k^3 = 0 \binom{n}{1} + 1 \binom{n}{2} + 6 \binom{n}{3} + 6 \binom{n}{4} = \lfrac{n(n-1)}{2} \Big( 1 + \lfrac{6(n-2)}{3} \big( 1 + \lfrac{n-3}{4} \big) \Big) = \Big( \lfrac{n(n-1)}{2} \Big)^2$.

Computation efficiency

This is far more efficient than the typical method in some textbooks (namely by taking summation on both sides of $(n+1)^3-n^3 = 3n^2+3n+1$ and telescoping) because the Newton series is easy to compute and manipulate. For sum of $p$-powers we only need $O(p^2)$ arithmetic operations to find the forward-differences and then $O(p^2)$ more to simplify the series form into a standard polynomial form. In contrast, the other method requires $O(p^3)$ arithmetic operations.

Answer 13

One technique is the Newton series, which permits termwise summation just as the Taylor series permits termwise integration. There is also the more powerful but less efficient method of indefinite summation, which is unsurprisingly the discrete version of integration by parts. $ \def\zz{\mathbb{Z}} \def\lbinom#1#2{{\large\binom{#1}{#2}}} $

---

Let $R$ be the right-shift operator defined as:

$R = ( \text{function $f$ on $\zz$} \mapsto ( \text{int $n$} \mapsto f(n+1) ) )$

Namely, for any function $f$ on $\zz$ and $n \in Z$, we have $R(f)(n) = f(n+1)$.

For convenience we also define the summation operator:

$Σ = ( \text{function $f$ on $\zz$} \mapsto ( \text{int $n$} \mapsto \sum_{k=0}^{n-1} f(k) ) )$

Then we have the important property that $ΔΣ(f) = f$ for any function $f$ on $\zz$, analogous to the fundamental theorem of calculus.

Now we can derive summation by parts from the product rule in the same manner as for the continuous version. Consider any functions $f,g$ on $\zz$:

$Δ(f·g)(n) = f(n+1) g(n+1) - f(n) g(n) = f(n+1) Δ(g)(n) - Δ(f)(n) g(n)$ for any $n \in \zz$.

And hence the discrete product rule (with the usual pointwise sum and product of functions):

$Δ(f·g) = R(f)·Δ(g) + Δ(f)·g$.

Now by substituting $f$ with $Σ(f)$ and taking summation on both sides we get summation by parts:

$Δ(Σ(f)·g) = R(Σ(f))·Δ(g) + f·g$.

$Σ(f·g) = Σ(f)·g - Σ(R(Σ(f))·Δ(g)) + c$ for some constant function $c$ on $\zz$.

---

For convenience we let "$Δ_n E(n)$" denote "$E(n+1)-E(n)$" and "$Σ_n E(n)$" denote "$\sum_{k=0}^{n-1} E(k)$".

We can therefore compute $\sum_{k=0}^{n-1} k^3$ as follows:

$\sum_{k=0}^{n-1} k^3 = Σ_n\Big(\lbinom{n}{0}·n^3\Big)$   [Recall this definition of binomial coefficients.]

$\ = \lbinom{n}{1}·n^3 - Σ_n\Big(\lbinom{n+1}{1}·(3n^2+3n+1)\Big) + c$, where $c$ is some constant (independent of $n$),

$\ = \lbinom{n}{1}·n^3 - \lbinom{n+1}{2}·(3n^2+3n+1) + Σ_n\Big(\lbinom{n+2}{2}·(6n+6)\Big) + c'$,

    where $c'$ is some (other) constant,

$\ = \lbinom{n}{1}·n^3 - \lbinom{n+1}{2}·(3n^2+3n+1) + \lbinom{n+2}{3}·(6n+6) - Σ_n\Big(\lbinom{n+3}{3}·6\Big) + c''$,

    where $c''$ is some (other) constant,

$\ = \lbinom{n}{1}·n^3 - \lbinom{n+1}{2}·(3n^2+3n+1) + \lbinom{n+2}{3}·(6n+6) - \lbinom{n+3}{4}·6 + c''$.

By substituting $n = 0$ we readily get $c'' = 0$ and hence simplifying gives the desired formula.

Answer 14

A bit late but here what I found: $$ \sum_{k=1}^n k^p = (n+1)\sum_{j=1}^n j^{p-1} - \sum_{k=1}^n \sum_{j=1}^{k} j^{p-1} $$ So $$ F(n,p)=\sum_{k=1}^n k^p $$ $$ F(n,p)= \begin{cases} (n+1)F(n,p-1)-\sum_{k=1}^n F(k,p-1), & \text{if $p>0$} \\ n,& \text{if $p=0$} \\ \end{cases} $$ I worked it out by rearranging the elements of the sum: $$ \sum_{k=1}^n k^p $$ $$ \sum_{k=1}^n k*k^{p-1} $$ $$ \sum_{k=1}^n \sum_{j=1}^k k^{p-1} $$ $$ \sum_{k=1}^n (\sum_{j=1}^n j^{p-1} - \sum_{j=1}^{k-1} j^{p-1}) $$ $$ n\sum_{j=1}^n j^{p-1} - \sum_{k=1}^n \sum_{j=1}^{k-1} j^{p-1} $$ $$ n\sum_{j=1}^n j^{p-1} - \sum_{k=1}^n (-k^{p-1} + \sum_{j=1}^{k} j^{p-1}) $$ $$ n\sum_{j=1}^n j^{p-1} + \sum_{k=1}^n k^{p-1} - \sum_{k=1}^n \sum_{j=1}^{k} j^{p-1} $$ $$ (n+1)\sum_{j=1}^n j^{p-1} - \sum_{k=1}^n \sum_{j=1}^{k} j^{p-1} $$ It's very inefficient but it shows a nice recursive relation and it's very easy to implement. An example in JavaScript:

const PowerSum = function(n, p) {
    if(p === 0)
        return n;
    let s = (n + 1) * PowerSum(n, p - 1);
    for(let k = 1; k <= n; k++)
        s -= PowerSum(k, p - 1);
    return s;
}