Sum of sum of $k$th power of first $n$ natural numbers.
$$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)$$
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.
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.