How to simplify this triple summation containing binomial coefficients? 
$$
  \large\sum_{i=0}^{n} \sum_{j=i}^{n} \sum_{k=j}^{n} 
    \binom{i+m-1}{m-1}\binom{j+m-1}{m-1}\binom{k+m-1}{m-1}
$$

How to solve it when this involve more than thousand summation ?
 A: Here's almost a solution.  Let the symmetric group $S_3$ act on the set of triples $\{(i,j,k)\mid 1\le i,j,k\le n\}$.  Let $f(i,j,k)$ denote the product of binomial coefficients being summed.  Note that $f$ is symmetric, so $f(i',j',k')=f(i,j,k)$ if $(i',j',k')$ and $(i,j,k)$ are in the same $S_3$ orbit.  Thus $f$ actually defines a function on the orbits.  The sum you are after is the sum of $f$ over all orbits.  This is an application of the cycle index formula.
Since the cycle index polynomial of $S_3$ is $\frac16(a_1^3+3a_2a_1+2a_3)$, we get your sum by evaluating this polynomial at $a_k=\sum_{i=0}^n\binom{i+m-1}{m-1}^k$.
The sum $a_1$ can be done by hand:
$$a_1=\frac{n+1}{m}\binom{m+n}{m-1}.$$
The other two are not so simple (since, for example, when $i=j$ you get the square of the corresponding binomial coefficient).  For fixed $m$, you can get a closed form (see below for $m=2$), but for general $m,n$ Mathematica gives
$$a_2=\frac{\Gamma (1-2 m)}{\Gamma (1-m)^2}-\binom{m+n}{m-1}^2 \,
   _3F_2(1,m+n+1,m+n+1;n+2,n+2;1)$$
and
$$\begin{multline}a_3=\frac{\left(1-\frac{m}{2}\right)_{-m} (m+1)_{-m}}{(1)_{-m}
   \left(\frac{m}{2}+1\right)_{-m}}-\\ \binom{m+n}{m-1}^3 \,
   _4F_3(1,m+n+1,m+n+1,m+n+1;n+2,n+2,n+2;1).
\end{multline}$$
Unfortunately, these expressions are indeterminate; I have not yet been able to coax a useful form out of Mathematica.
In the special case $m=2$, here are the closed forms:
$$\begin{align}
a_1=&\frac{1}{2} (n+1) (n+2) \\
a_2=&\frac{1}{6} (n+1) (n+2) (2 n+3) \\
a_3=&\frac{1}{4} (n+1)^2 (n+2)^2 \\
\frac16(a_1^3+3a_2a_1+2a_3)=&\frac{1}{48} (n+1)^2 (n+2)^2 (n+3) (n+4)
\end{align}$$
This gives the following values for small $n$ (which are easily verified to agree with the given summation):
\begin{array}{cc}
 n & \\
\hline
 0 & 1 \\
 1 & 15 \\
 2 & 90 \\
 3 & 350 \\
 4 & 1050 \\
 5 & 2646 \\
\end{array}
A: Let $m, n, s$ stand for any positive integers. In following discussion, we will assume $m$ is given.
Let $[n] = \{\; 0, 1, \ldots, n \;\}$  and define
$$a(n) = \binom{m-1+n}{m-1}
\quad\text{ and }\quad
f_s(n) = \sum_{k=1}^n a(n)^s
$$
The sum $\mathcal{S}$ we want  can be rewritten as
$$\mathcal{S} = \sum_{0 \le i\le j \le k \le n} a(i)a(j)a(k)$$
Start with any triple $(i,j,k)$ appear in above sum, if we reorder them, they will generate every element of $[n]^3$ at least once. In general,


*

*if $i \ne j \ne k$, it will generate 6 distinct triples.

*If two indices coincides, say $i = j \ne k$ or $i \ne j = k$, it will generate three distinct triples.

*If $i = j = k$, there is only one distinct triples.


As a consequence of this, we find
$$\begin{align}
6 \mathcal{S} &= f_1(n)^3 + (6-3) f_1(n)f_2(n) + (6-3-1)f_3(n)\\
\iff
\mathcal{S} &= \frac16\big( f_1(n)^3 + 3 f_1(n)f_2(n) + 2f_3(n)\big)
\end{align}
$$
Our task reduces to the computation of $f_1(n)$, $f_2(n)$ and $f_3(n)$.
For $f_1(n)$, it is easy. Since the OGF for $a(n)$ is
$$\sum_{n=0}^\infty a(n) z^n = \sum_{n=0}^\infty \binom{m-1+n}{m-1} z^n = \frac{1}{(1-z)^m}$$
The OCF for $f_1(n)$ will be equal to
$$\sum_{n=0}^\infty f_1(n) z^n = \frac{1}{(1-z)(1-z)^{m}} =
\frac{1}{(1-z)^{m+1}}$$
This leads to $\displaystyle\;f_1(n) = \binom{m+n}{m}\;$.
For other $f_s(n)$ with $s > 1$, we will use the fact $a(n)$, as a function of $n$, is a polynomial in $n$ with degree $m-1$. This means $a(n)^s$ itself is a polynomial in $n$ with degree $s(m-1)$ and hence satisfies a linear recurrence relation of degree $s(m-1)+1$. Translate this into OCF, we find
$$\sum_{n=0}^\infty a(n)^s z^n = \frac{B_s(z)}{(1-z)^{s(m-1)+1}}$$
for some polynomial $B_s(z)$ with degree $\le s(m-1)$. We can backout the coefficients of $B_s(z)$ by comparing the two sides of following relations
$$B_s(z) = (1-z)^{s(m-1)+1} \left(\sum_{n=0}^\infty a(n)^s z^s\right)$$
The end result is
$$B_s(z) = \sum_{n=0}^{s(m-1)} b_s(n) z^n
\;\text{ with }\;
b_s(n) = \sum_{k=0}^n (-1)^k \binom{s(m-1)+1}{k} a(n-k)^s
$$
As a result, the OGF for $f_s(n)$ is
$$ \frac{\sum\limits_{n=0}^{s(m-1)} b_s(n) z^n}{(1-z)^{s(m-1)+2}}
\quad\implies\quad
f_s(n) = \sum_{k=0}^{s(m-1)}b_s(k)\binom{s(m-1)+1+n-k}{s(m-1)+1}
$$
With a little bit of algebra, one can simplify $b_2(n)$ and $b_3(n)$ as
$$\left\{\begin{align}
b_2(n) &= \binom{m-1}{n}^2\\
b_3(n) &= \sum_{k=\max(0,n-m+1)}^{\min(n,m-1)} \binom{m-1}{k}^2\binom{m-1+k}{n}\binom{2m-2-k}{2m-2-n}
\end{align}\right.$$
and in above sum, $b_2(n)$ and $b_3(n)$, terminate earlier at $m-1$ and $2(m-1)$ respectively.
As a final summary, we have
$$
\mathcal{S} = \frac16\big( f_1^3 + 3 f_1f_2 + 2f_3\big)
\quad\text{ with }\quad
\left\{\begin{align}
f_1 &= \binom{m+n}{m}\\
f_2 &= \sum_{k=0}^{m-1} b_2(k)\binom{2m-1+n-k}{2m-1}\\
f_3 &= \sum_{k=0}^{2m-2}b_3(k)\binom{3m-2+n-k}{3m-2}\\
\end{align}\right.
$$
Using this, one can compute $S$ for small $m$ as a polynomial in $n$.


*

*$m = 1$, $S = \binom{n+3}{3}$.

*$m = 2$, $S = \binom{n+4}{4}\binom{n+2}{2}$.

*$m = 3$, $S = \binom{n+5}{5}\frac{35n^4+294n^3+844n^2+969n+378}{378}$.

*$m = 4$, $S = \binom{n+6}{6}\binom{n+4}{4}\frac{35n^2+205n+168}{168}$.

A: A second answer from me, the first double/triple sum free version,
verified up to a certain degree.
display2d : false;

for n : 0 thru 50 do block(
        print(n),
for m : 0 thru 50 do block(
        mylhs : sum(sum(sum(
                    1
                    *binomial(i+m,m)
                    *binomial(j+m,m)
                    *binomial(k+m,m)
        ,k,j,n),j,i,n),i,0,n),
        myrhs : 0,
        myrhs : myrhs + (
                    1
                    *binomial(n+m+1,m+1)
                    *binomial(n+m+1,m+1)
                    *binomial(n+m+1,m+1)
        ),
        myrhs : myrhs - sum(
                    1
                    *binomial(i+m,m)
                    *binomial(i+m,m+1)
                    *binomial(n+m+1,m+1)
        ,i,0,n),
        myrhs : myrhs - sum(
                    1
                    *binomial(i+m,m)
                    *binomial(i+m,m+1)
                    *binomial(i+m+1,m+1)
        ,i,0,n),
        res : mylhs/myrhs,
        if (res # 1) then block(
                print(n,m,"|",res)
        )
));

