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

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.

Solution 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.

Solution 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)!}$.