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