# Sum of sum of $k$th power of first $n$ natural numbers.

I was working on a problem which involves computation of $k$-th power of first $n$ natural numbers.

Say $f(n) = 1^k+2^k+3^k+\cdots+n^k$ we can compute $f(n)$ by using Faulhaber's Triangle also by fast computation of Bernoulli numbers.

Here we have to compute $G(n)= f(1)+f(2)+f(3)+\cdots+f(n)$

I did some work on this to simplify but unable to come up with any easy equation.

Can someone help me how to compute $G(n)$ efficiently , I can efficiently compute $f(n)$.

Here $n<123456789$ and $k<321$.

Thanks.

• By your problem, you have $f(i)=1^i+2^i+\cdots+n^i$. Then- $G(n)$ $=$ $\displaystyle \sum_{i=1}^n f(i)$ $=$ $\displaystyle \sum_{i=1}^n 1^i$ $+$ $\displaystyle \sum_{i=1}^n 2^i$ $\cdots$ $+$ $\displaystyle \sum_{i=1}^n n^i$ $=$ $n$ $+$ $\displaystyle \sum_{i=2}^n\dfrac {i(i^n-1)}{i-1}$. – William Hilbert May 20 '14 at 16:57
• There are formulas using Bernoulli numbers. – vonbrand May 20 '14 at 16:58
• Are those formulas really helpful in solving this problem? – William Hilbert May 20 '14 at 17:00
• @William Hilbert its not like that given n and k we have to compute f(n)=1^k+2^k+3^k...n^k which is completely different than what you are saying. – Luckymaster May 20 '14 at 17:01
• Then perhaps I misunderstood the problem. Can you tell me what the problem says? – William Hilbert May 20 '14 at 17:02

$$G_k(n)= \sum_{i=1}^n f_k(i)= \sum_{i=1}^n \sum_{j=1}^i j^k = \sum_{j=1}^n (n+1-j) j^k=(n+1)\sum_{j=1}^n j^k - \left(\sum_{j=1}^n j^{k+1} \right)\\ =(n+1)f_k(n)-f_{k+1}(n)$$

Since you said you can efficiently calculate

$$f_{k}(n) =\sum_{j=1}^n j^k$$ you can also calculate efficiently $$(n+1)f_k(n)-f_{k+1}(n)=g_k(n)$$

• Ah yeah! It was fairly trivial. :) – Piotr Miś May 20 '14 at 18:02
• @ N.S thank you sir. this is what I was looking for and thanks all for their great answers. – Luckymaster May 21 '14 at 4:24

One more type of answer is to use the Faulhaber-matrix (with a little tweak) and iterate.
Let's denote one vector of consecutive powers of $x$ as $V(x)=[1,x,x^2,x^3,...]^T$ either up to required dimension (see below) or simply assume as of infinite size. Then the dotproduct with the (tweaked) Faulhabermatrix with that vector gives the vector of sums-of-like powers.

The top-left of the Faulhabermatrix, with one empty column prepended:
$$G= \small \begin{bmatrix} . & 1 & . & . & . & . & . & . \\ . & 1/2 & 1/2 & . & . & . & . & . \\ . & 1/6 & 1/2 & 1/3 & . & . & . & . \\ . & . & 1/4 & 1/2 & 1/4 & . & . & . \\ . & -1/30 & . & 1/3 & 1/2 & 1/5 & . & . \\ . & . & -1/12 & . & 5/12 & 1/2 & 1/6 & . \\ . & 1/42 & . & -1/6 & . & 1/2 & 1/2 & 1/7 \\ . & . & 1/12 & . & -7/24 & . & 7/12 & 1/2 \end{bmatrix} \tag 1$$ Then $$G \cdot V(x)=S(x) \tag 2$$ where due to Faulhabers'analysis $$S(x)=V(1)+V(2)+V(3)+...+V(x) \\ \qquad \qquad \quad=[x, 1+2+...x,1^2+2^2+...+x^2,...]^T \tag 3$$
Iterating this gives a fairly intuitive method do arrive at the "sum-of-sums" $SS(x)$ $$\begin{array} {} G \cdot S(x) &= SS(x) \\ &= G\cdot V(1)+ G\cdot V(2)+...+G \cdot V(x) \\ &= V(1)+ \left(V(1)+V(2)\right)+...+\left(V(1)+V(2)+...+V(x)\right) \\ &= [1+1+1+...+1, \\ & \qquad 1+(1+2)+(1+2+3)+...(1+2+3+...+x) , \\ & \qquad 1+(1+2^2)+(1+2^2+3^2)+...(1+2^2+3^2+...+x^2) , \\ & \qquad ... , \\ & \qquad1+(1+2^k)+(1+2^k+3^k)+...(1+2^k+3^k+...+x^k), \\ & \qquad ...]^T \\\end{array} \tag 4$$ This can be rewritten as $$G^2 \cdot V(x) = SS(x) \tag 5$$ and thus that square of the Faulhaber-matrix is of interest. We find the top-left of $G^2$ : $$G^2 = \small \begin{bmatrix} . & 1/2 & 1/2 & . & . & . & . & . \\ . & 1/3 & 1/2 & 1/6 & . & . & . & . \\ . & 1/6 & 5/12 & 1/3 & 1/12 & . & . & . \\ . & 1/30 & 1/4 & 5/12 & 1/4 & 1/20 & . & . \\ . & -1/30 & 1/20 & 1/3 & 5/12 & 1/5 & 1/30 & . \\ . & -1/42 & -1/12 & 1/12 & 5/12 & 5/12 & 1/6 & 1/42 \\ . & 1/42 & -5/84 & -1/6 & 1/8 & 1/2 & 5/12 & 1/7 \\ . & 1/30 & 1/12 & -5/36 & -7/24 & 7/40 & 7/12 & 5/12 \end{bmatrix} \tag 6$$ and indeed, the dot-product with $V(3)$ looks like $$G^2 \cdot V(3) = SS(3) = \small \begin{bmatrix} 6 \\ 10 \\ 20 \\ 46 \\ 116 \\ 310 \\ 860 \\ 2446 \end{bmatrix} \tag 7$$ where in the $k$'th row is the required sum with summands to the $k$'th power.

Using dimension $8$ only: $$G^2 \cdot V(123456789) = \Tiny \begin{bmatrix} 7620789436823655 \\ 313612736252315226397035 \\ 19358810860413733959550603308825 \\ 1433985988226290537213517145779502623643 \\ 118023538007597075768233960814021641290827487705 \\ 10407719390614949736958778538255639230572384078931269755 \\ 963677720389558291637495115742280128808534343727190032990523625 \\ 92534211741854093396678577789655310126697035147135144478399613856675643 \end{bmatrix}$$ Using dimension $n=322$ (so that the highest exponent is $k=321$ in row $322$)

n=322
result = Gp^2 \cdot V(123456789)

k  | ss_k(123456789)
---!----------------------
0    7.62078943682E15
1    3.13612736252E23
2    1.93588108604E31
3    1.43398598823E39
4    1.18023538008E47
...
317  1.53837082910E2576
318  1.88735309863E2584
319  2.31554800916E2592
320  2.84094533467E2600
321  3.48562264082E2608

## Generalization

That use of the $G$-matrix allows immediate generalization to further iterates by further powers-of-$G$. Even more, we have the simple option to let our sums start at another value than $[1,1,1^2,1^3,...]^T$ for instance $$G^2 \cdot \left( V(5)-V(2) \right) = SS(5)-SS(2) = S(5)+S(4)+S(3) = \small \begin{bmatrix} 12 \\ 31 \\ 99 \\ 361 \\ 1431 \\ 6001 \\ 26199 \\ 117841 \end{bmatrix} \tag 8$$

## Appendix:

If you want to exercise with this it might be easier to create the matrix $G$ from the Faulhaber-matrix $F$ and $F$ itself as inversion or the reduced Pascalmatrix: $$P^* = \small \left[ \begin{array} {r} 1 & . & . & . & . & . & . & . \\ -1 & 2 & . & . & . & . & . & . \\ 1 & -3 & 3 & . & . & . & . & . \\ -1 & 4 & -6 & 4 & . & . & . & . \\ 1 & -5 & 10 & -10 & 5 & . & . & . \\ -1 & 6 & -15 & 20 & -15 & 6 & . & . \\ 1 & -7 & 21 & -35 & 35 & -21 & 7 & . \\ -1 & 8 & -28 & 56 & -70 & 56 & -28 & 8 \end{array} \right]$$ which is easier to create programmatically and then $$F = (P^*)^{-1} \tag 9$$ and $G$ is $F$ with one empty column prepended.

You might find it useful to visualize the sum: $$\begin{matrix} 1^k \\ 1^k & 2^k \\ 1^k & 2^k & 3^k \\ 1^k & 2^k & 3^k & 4^k \\ \vdots&\vdots&\vdots&\vdots&\ddots \\ 1^k & 2^k & 3^k & 4^k & \dots & n^k \\ \end{matrix}$$

$$G_k(n) = n\cdot1^k + (n-1)\cdot2^k + (n-2)\cdot3^k + \dots + 2\cdot(n-1)^k + 1\cdot n^k = \sum_{s=1}^{n} (n-s+1)s^k$$

This will take $O(n\lg k)$ multiplications to compute directly. I'm trying to find a better way.

EDIT: I was on the right track. For glorious ending see answer posted by N. S.

You can compute the sum of first $n$ $k$-th powers in $O(k^2)$ We define $$F_k(n) = \sum_{i=1}^{n} i^k$$ Forse example, $F_0(n) = n$ and $F_1(n) = \frac{n(n+1)}{2}$. Now, we can prove that $$F_k(n) = (n+1)^{k+1} -1 - \sum_{i=2}^{k+1} \binom{k+1}{i} \cdot F_{k-1+i}(n)$$ You simply compute the first $k$ rows of Pascal Triangle, and then for every $j \le k$ you compute $F_j(n)$