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.

Solution 1:

$$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)$$

Solution 2:

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.

Your example:

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.

Solution 3:

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.