Calculate the sum of first $n$ natural numbers taken $k$ at a time
Solution 1:
Overview: This answer consists of three parts
At first we develop summation formulae for slightly different $k$-fold sums for $k=1,2,3$. We take indices starting from $0$ instead of $1$, which makes calculation somewhat more convenient.
We obtain for general $k\geq 1$ a nice expression of the $k$-fold sum \begin{align*} \sum_{i_1=0}^{n}\sum_{i_2=0}^{n-i_1}&\ldots\sum_{i_k=0}^{n-i_1-i_2-\ldots-i_{k-1}}i_1(i_1+i_2)\cdots(i_1+i_2+\ldots+i_k)\\ \end{align*} as sum of binomial coefficients .
We look at the difference of the $k$-fold sums when starting from $0$ and when starting from $1$.
Hint: If you are not familiar with generating functions you might want to look at the backstage info at the end of this answer.
Part 1: Summation formula for $k=1,2,3$ and start index $0$
The following is valid
\begin{align*} \sum_{i=0}^ni&=\binom{n+1}{2}=\frac{1}{2}n(n+1)\tag{1}\\ \sum_{i=0}^n\sum_{j=0}^{n-i}i(i+j)&=\binom{n+3}{4}+2\binom{n+2}{4}\\ &=\frac{1}{24}n(n+1)(n+2)(3n+1)\tag{2}\\ \sum_{i=0}^n\sum_{j=0}^{n-i}\sum_{k=0}^{n-i-j}i(i+j)(i+j+k)&=\binom{n+5}{6}+8\binom{n+4}{6}+6\binom{n+3}{6}\\ &=\frac{1}{48}n^2(n+1)^2(n+2)(n+3)\tag{3}\\ \end{align*}
Note, the representation of binomial coefficients looks promising, cause a pattern for $k>3$ can easily be derived. We will see in Part 2, that the multiplicities of the binomial coefficients also follow a well-known pattern.
In order to show these formulae we need a repertoire of $\left(x\frac{d}{dx}\right)^k\frac{1}{1-x}$ for $k=1,\ldots,3$.
We obtain \begin{align*} \left(x\frac{x}{dx}\right)\frac{1}{1-x}&=\frac{x}{(1-x)^2} =\sum_{n=0}^{\infty}nx^n\\ \left(x\frac{x}{dx}\right)^2\frac{1}{1-x}&=\frac{x(1+x)}{(1-x)^3} =\sum_{n=0}^{\infty}n^2x^n\\ \left(x\frac{x}{dx}\right)^3\frac{1}{1-x}&=\frac{x(1+4x+x^2)}{(1-x)^4} =\sum_{n=0}^{\infty}n^3x^n\\ \end{align*}
It's also convenient to use the coefficient of operator $[x^n]$ to denote the coefficient of $x^n$ in a generating series.
We obtain for $k=1$ \begin{align*} \sum_{i=0}^{n}i&=\sum_{i=0}^{n}i\cdot1\\ &=\sum_{i=0}^n[x^i]\left(x\frac{x}{dx}\right)\frac{1}{1-x}[x^{n-i}]\frac{1}{1-x}\tag{4}\\ &=\sum_{i=0}^n[x^i]\frac{x}{(1-x)^2}[x^{n-i}]\frac{1}{1-x}\\ &=[x^n]\frac{x}{(1-x)^3}\tag{5}\\ &=[x^n]x\sum_{k=0}^{\infty}\binom{-3}{k}(-x)^k\\ &=[x^{n-1}]\sum_{k=0}^{\infty}\binom{k+2}{2}x^k\tag{6}\\ &=\binom{n+1}{2}\\ &=\frac{1}{2}n(n+1) \end{align*} and (1) follows.
Comment:
In (4) we use $$i=[x^i]\sum_{k=0}^{\infty}kx^x=[x^i]\left(x\frac{x}{dx}\right)\frac{1}{1-x}$$ and $$1=[x^{n-i}]\sum_{k=0}^{\infty}x^k=[x^{n-i}]\frac{1}{1-x}$$
In (5) we use $[x^n]A(x)B(x)=\sum_{i=0}^{n}\left([x^i]A(x)\right)\left([x^{n-i}]B(x)\right)$
In (6) we use the formula $\binom{-n}{k}=\binom{n+k-1}{k}(-1)^k=\binom{n+k-1}{n-1}(-1)^k$
We obtain for $k=2$ \begin{align*} \sum_{i=0}^{n}&\sum_{j=0}^{n-i}i(i+j)\\ &=\sum_{i=0}^ni^2\sum_{j=0}^{n-i}1+\sum_{i=0}^ni\sum_{j=0}^{n-i}j\\ &= \sum_{i=0}^n[x^i]\left(x\frac{x}{dx}\right)^2\frac{1}{1-x}[x^{n-i}]\frac{1}{(1-x)^2}\\ &\qquad+\sum_{i=0}^n[x^i]\left(x\frac{x}{dx}\right)\frac{1}{1-x}[x^{n-i}]\frac{1}{1-x}\left(x\frac{x}{dx}\right)\frac{1}{1-x}\\ &= \sum_{i=0}^n[x^i]\frac{x(1+x)}{(1-x)^3}[x^{n-i}]\frac{1}{(1-x)^2}+\sum_{i=0}^n[x^i]\frac{x}{(1-x)^2}[x^{n-i}]\frac{x}{(1-x)^3}\\ &=[x^n]\frac{x+2x^2}{(1-x)^5}\\ &=\left([x^{n-1}]+2[x^{n-2}]\right)\sum_{k=0}^{\infty}\binom{-5}{k}(-x)^k\\ &=\left([x^{n-1}]+2[x^{n-2}]\right)\sum_{k=0}^{\infty}\binom{k+4}{4}x^k\\ &=\binom{n+3}{4}+2\binom{n+2}{4}\\ &=\frac{1}{24}n(n+1)(n+2)(3n+1) \end{align*} and (2) follows.
$$ $$
We obtain for $k=3$ \begin{align*} \sum_{i=0}^{n}&\sum_{j=0}^{n-i}\sum_{k=0}^{n-i-j}i(i+j)(i+j+k)\\ &=\sum_{i=0}^{n}\sum_{j=0}^{n-i}\sum_{k=0}^{n-i-j}(i^3+2i^2j+i^2k+ij^2+ijk)\\ &=\sum_{i=0}^{n}\sum_{j=0}^{n-i}\sum_{k=0}^{n-i-j}(i^3+4i^2j+ijk)\tag{7}\\ &=\sum_{i=0}^{n}i^3\sum_{j=0}^{n-i}1\sum_{k=0}^{n-i-j}1 +4\sum_{i=0}^{n}i^2\sum_{j=0}^{n-i}j\sum_{k=0}^{n-i-j}1\tag{8}\\ &\qquad+\sum_{i=0}^{n}i\sum_{j=0}^{n-i}j\sum_{k=0}^{n-i-j}k\\ \end{align*}
Note, in (7) we use the symmetry \begin{align*} \sum_{i=0}^{n}\sum_{j=0}^{n-i}\sum_{k=0}^{n-i-j}i^2j =\sum_{i=0}^{n}\sum_{j=0}^{n-i}\sum_{k=0}^{n-i-j}i^2k =\sum_{i=0}^{n}\sum_{j=0}^{n-i}\sum_{k=0}^{n-i-j}ij^2 \end{align*}
The calculation of the three sums in (8) is straight forward and can be done similarly to $k=1,2$.
We obtain
\begin{align*} \sum_{i=0}^{n}i^3\sum_{j=0}^{n-i}1\sum_{k=0}^{n-i-j}1&=[x^n]\frac{x(1+4x+x^2)}{(1-x)^7} =\binom{n+3}{6}+4\binom{n+4}{6}+\binom{n+5}{6}\\ \sum_{i=0}^{n}i^2\sum_{j=0}^{n-i}j\sum_{k=0}^{n-i-j}1&=[x^n]\frac{x^2(1+x)}{(1-x)^7} =\binom{n+4}{6}+\binom{n+3}{6}\\ \sum_{i=0}^{n}i\sum_{j=0}^{n-i}j\sum_{k=0}^{n-i-j}k&=[x^n]\frac{x^3}{(1-x)^7} =\binom{n+3}{6}\\ \end{align*}
Combining the three sums according to (7) results in
\begin{align*} \sum_{i=0}^{n}\sum_{j=0}^{n-i}\sum_{k=0}^{n-i-j}(i^3+4i^2j+ijk) &=6\binom{n+3}{6}+8\binom{n+4}{6}+\binom{n+5}{6}\\ &=\frac{1}{48}n^2(n+1)^2(n+2)(n+3) \end{align*} and the claim (3) follows.
$$ $$
Part 2: Summation formula for all $k\geq 1$ and startindex $0$.
In the following we don't give a proof for general $k$, but we provide some aspects which give strong evidence for the correctness of the claim.
When looking at
\begin{align*} \sum_{i=0}^ni&=\color{blue}{1}\binom{n+1}{2}\\ \sum_{i=0}^n\sum_{j=0}^{n-i}i(i+j)&=\color{blue}{1}\binom{n+3}{4}+\color{blue}{2}\binom{n+2}{4}\\ \sum_{i=0}^n\sum_{j=0}^{n-i}\sum_{k=0}^{n-i-j}i(i+j)(i+j+k)&=\color{blue}{1}\binom{n+5}{6}+\color{blue}{8}\binom{n+4}{6}+\color{blue}{6}\binom{n+3}{6}\\ \end{align*} the shape of the binomial coefficients can be easily generalized. The coefficients \begin{align*} \color{blue}{1};\quad \color{blue}{1},\color{blue}{2};\quad \color{blue}{1},\color{blue}{8},\color{blue}{6} \end{align*} are part of the OEIS sequence A008517 and give the values of the Second order Eulerian Triangle $T(k,l),1\leq l\leq k$.
The values $T(4,l), 1\leq l\leq 4$ are $\color{blue}{1},\color{blue}{22},\color{blue}{58},\color{blue}{24}$ and indeed, it is easy to verify that
\begin{align*} \sum_{i=0}^n&\sum_{j=0}^{n-i}\sum_{k=0}^{n-i-j}\sum_{l=0}^{n-i-j-l}i(i+j)(i+j+k)(i+j+k+l)\\ &= \color{blue}{1}\binom{n+7}{8}+\color{blue}{22}\binom{n+6}{8}+\color{blue}{58}\binom{n+5}{8}+\color{blue}{24}\binom{n+4}{8} \end{align*}
We are now in the position to
Claim: The following formula of the $k$-fold sum is valid \begin{align*} \sum_{i_1=0}^{n}\sum_{i_2=0}^{n-i_1}&\ldots\sum_{i_k=0}^{n-i_1-i_2-\ldots-i_{k-1}}i_1(i_1+i_2)\cdots(i_1+i_2+\ldots+i_k)\\ &=\sum_{l=1}^{k}\color{blue}{T(k,l)}\binom{n+2k-l}{2k}\tag{9} \end{align*} with $T(k,l)$ the numbers of the second order Eulerian Triangle (A008517).
$$ $$
Part 3: Summation formula for all $k\geq 1$ and startindex $1$.
Finally some aspects about the summation formula with indices starting with $1$.
When looking at the differences of \begin{align*} \sum_{i=0}^ni&-\sum_{i=1}^ni=0\\ \sum_{i=0}^n\sum_{j=0}^{n-i}i(i+j)&-\sum_{i=1}^n\sum_{j=1}^{n-i}i(i+j)=\frac{1}{6}n(n+1)(2n+1)\\ \sum_{i=0}^n\sum_{j=0}^{n-i}\sum_{k=0}^{n-i-j}i(i+j)(i+j+k)&- \sum_{i=1}^n\sum_{j=1}^{n-i}\sum_{k=1}^{n-i-j}i(i+j)(i+j+k)\\ &\qquad=\frac{1}{12}(n+1)^2(n+2)^2(2n+3)\\ \end{align*}
we can find them as A000330 ($k=2)$ and as A108674 ($k=3)$ in the OEIS database.
Some further elaboration could be to look for a difference formula for general $k$ and combine it with the result of (9) to obtain a nice expression for OPs $k$-fold sum with indices starting from $1$.
Backstage info:
You might want to skip this section, if generating functions are already known.
When looking at the $2$-fold sum \begin{align*} \sum_{i=1}^n\sum_{j=1}^{n-i}i(i+j)=\sum_{i=1}^n\left(i^2\sum_{j=1}^{n-i}1\right)+\sum_{i=1}^n\left(i\sum_{j=1}^{n-i}j\right) \end{align*} the outer sum has the shape of a Cauchy product \begin{align*} \sum_{i=0}^{n}a_ib_{n-i} \end{align*} with index $i$ starting at $1$ instead of $0$.
Since we want to use generating functions $A(x)=\sum_{n= 0}^{\infty}a_nx^n$ to derive the summation formulae, and the product of generating function
\begin{align*} A(x)B(x)=\sum_{k= 0}^{\infty}a_kx^k\sum_{l= 0}^{\infty}b_lx^l = \sum_{n= 0}^{\infty}\left(\sum_{k=0}^{n-k}a_kb_{n-k}\right)x^n \end{align*} gives Cauchy products, we consider $k$-fold sums with index starting from $0$ instead.
This is interesting by its own and later we can look at the difference to $k$-fold sums with indices starting from $1$.
We can successively apply the $\left(x\frac{d}{dx}\right)$-operator to a generating function $A(x)$ to obtain
\begin{align*} \left(x\frac{d}{dx}\right)A(x)&=\sum_{n=0}^{\infty}na_nx^n\\ \left(x\frac{d}{dx}\right)^2A(x)&=\sum_{n=0}^{\infty}n^2a_nx^n \end{align*}
Multiplication of $A(x)$ with $\frac{1}{1-x}$ results in summing up the coefficients $a_n$ \begin{array}{crl} (a_n)_{n\geq 0}\qquad &\qquad A(x)=&\sum_{n=0}^{\infty}a_nx^n\\ \left(\sum_{k=0}^{n}a_k\right)_{n\geq 0}\qquad&\qquad\frac{1}{1-x}A(x)=&\sum_{n=0}^{\infty}\left(\sum_{k=0}^{n}a_k\right)x^n \end{array}