A Vandermonde Identity for Stirling Numbers? I'm facing the problem of trying to express a quantity in the simplest possible way (it is, using the least possible number of sum symbols).
$$ \sum_{j=0}^{n} \sum_{\ell=0}^m \frac{1}{j!}\binom{b+j}{j} {j+1 \brack {\ell+1}} {b+2 \brack {m-\ell+1}}$$
Of course, this can be easily written as a convolution between two polynomials (which happen to be more or less simple). I'm pretty sure that approach will not work (at most, one can write the above expression as "the coefficient of $x^m$ in this product [...]", but that is not useful to my purpose).
However, if one explores this sum a little bit, it pretty soon come up the fact that it could be truly useful to, for example, be able to compute this: 
$$\sum_{\ell=0}^m {j+1\brack{\ell+1}}{b+2 \brack {m-\ell+1}}$$
(which resembles a lot Vandermonde's Identity, but with Stirling numbers instead of binomial coefficients).
I looked up on a couple of books (Concrete Mathematics of Graham-Knuth-Patashnik, and others), and I couldn't find any references pointing to such an identity. Does anybody know something like that? (Perhaps involving other weird numbers as Eulerian or double Eulerian or that kind of stuff?)
Nevertheless, any kind of help simplifying the first double sum would be really appreciated.
 A: We show  the following identity     is   valid  for non-negative integers  $a,b,M$:
\begin{align*}     
\sum_{k=0}^M       {a\brack         k}{b\brack  M-k}=\sum_{n=0}^{\min{\{a,b\}}}{a+b-n\brack  M}\frac{(-1)^na!b!}{n!(a-n)!(b-n)!}\tag{1}
\end{align*}
The right-hand sum is regrettably not a simplification. But we have at least one factor as Stirling number of the   first kind instead of two (at the cost  of some other factors) and we have an  upper index  containing $a+b$ which comes somewhat close to Vandermonde's identity
\begin{align*}
\sum_{k=0}^m\binom{a}{k}\binom{b}{M-k}=\binom{a+b}{M}
\end{align*}
We  show (1) by recalling the generating function:
\begin{align*}
(1-z)^{-u}=\sum_{n=0}^\infty\sum_{k=0}^n{n\brack k}u^k\frac{z^n}{n!}\tag{2}
\end{align*}

At first we derive the left-hand side of (1). We use the coefficient of operator $[z^n]$ to denote the coefficient of $z^n$ and  obtain from (2) 
  \begin{align*}
\color{blue}{a!b!}&\color{blue}{[u^Mz^at^b](1-z)^{-u}(1-t)^{-u}}\tag{3}\\
&=a!b![u^M]\left([z^a]\sum_{n=0}^\infty\sum_{k=0}^n{n\brack k}u^k\frac{z^n}{n!}\right)\left([t^b]\sum_{n=0}^\infty\sum_{k=0}^n{n\brack k}u^k\frac{t^n}{n!}\right)\\
&=a!b![u^M]\left(\sum_{k=0}^a{a\brack k}u^k\frac{1}{a!}\right)\left(\sum_{l=0}^b{b\brack l}u^l\frac{1}{b!}\right)\\
&=[u^M]\sum_{q=0}^{a+b}\sum_{{k+l=q}\atop{k,l\geq 0}}{a\brack k}{b\brack l}u^q\\
&=[u^M]\sum_{q=0}^{a+b}\sum_{k=0}^q{a\brack k}{b\brack q-k}u^q\\
&\,\,\color{blue}{=\sum_{k=0}^M{a\brack k}{b\brack M-k}}\tag{4}
\end{align*}

We take (3) and calculate the coefficient in a somewhat different way.

We obtain
  \begin{align*}
\color{blue}{a!b!}&\color{blue}{[u^Mz^at^b]((1-z)(1-t))^{-u}}\\
&=a!b![u^Mz^at^b](1-(z(1-t)+t))^{-u}\\
&=a!b![u^Mz^at^b]\left(\sum_{n=0}^\infty\sum_{k=0}^n{n\brack k}u^k\frac{(z(1-t)+t)^n}{n!}\right)\tag{5}\\
&=a!b![z^at^b]\sum_{n=0}^\infty{n\brack M}\frac{(z(1-t)+t)^n}{n!}\tag{6}\\
&=a!b![z^at^b]\sum_{n=0}^\infty{n\brack M}\frac{1}{n!}\sum_{j=0}^n\binom{n}{j}z^j(1-t)^jt^{n-j}\tag{7}\\
&=a!b![t^b]\sum_{n=0}^\infty{n\brack M}\frac{1}{n!}\binom{n}{a}(1-t)^at^{n-a}\tag{8}\\
&=a!b!\sum_{n=0}^{\min\{a,b\}}{n\brack M}\frac{1}{n!}\binom{n}{a}\binom{a}{a+b-n}(-1)^{a+b-n}\tag{9}\\
&=a!b!\sum_{n=0}^{\min\{a,b\}}{n\brack M}\frac{(-1)^{a+b-n}}{(n-a)!(n-b)!(a+b-n)!}\tag{10}\\
&\,\,\color{blue}{=a!b!\sum_{n=0}^{\min{\{a,b\}}}{a+b-n\brack M}\frac{(-1)^n}{n!(a-n)!(b-n)!}}\tag{11}
\end{align*}
  and the claim (1) follows.

Comment:


*

*In (5) we expand the series according to (2).

*In (6) we select the coefficient of $u^M$.

*In (7) we expand the binomial.

*In (8) we select the coefficient of $z^a$.

*In (9) we select the coefficient of $t^b$.

*In (10) we do some cancellations.

*In (11) we change the order of summation $n\to a+b-n$.
A: As usual check carefully (although there is a code snippet to illustrate).
I seem to have a form of mathematical dyslexia.
Answer: $\left[x^{M}\right]\frac{\Gamma\left(x+a\right)\Gamma\left(x+b\right)}{\Gamma\left(x\right)\Gamma\left(x+a+b\right)}\left(x\right)_{\left(a+b\right)}$
Actually the expression gives the answer for all $M$. Some code is included for that. I can expand it to include M and the binomial summation on a if needed. 
The process can probably be specialized to just compute one answer for a particular $M$.  Let me know if you need it.
First some standard facts:
Definition 1. The Pochhammer symbol.
$$\begin{align*}\left(x\right)_{n}=x\cdot\left(x+1\right)\ldots\left(x+n-1\right)
\end{align*}$$ 
Definition 2. Unsigned Stirling Number
$$\left[\begin{array}{c}
n\\
l
\end{array}\right]\equiv\left[x^{l}\right]{\displaystyle \prod_{k=0}^{n-1}\left(x+k\right)=\left[x^{l}\right]x\cdot\left(x+1\right)\ldots\left(x+n-1\right)}=\left[x^{l}\right]\left(x\right)_{n}$$
$$=\left[x^{l}\right]{\displaystyle \sum_{j=0}^{n}}\left[\begin{array}{c}
n\\
j
\end{array}\right]x^{j}$$
We can partition  
$$\left(x\right)_{\left(a+b\right)}\rightarrow\left(x\right)_{a}\cdot\left(x+a\right)_{b} \tag{1}\label{1}$$
The formula for all convolutions:
$$(x)_{a}\cdot\left(x\right)_{b} \tag{2}\label{2}$$
The problem can be phrased as the convolution : 
$${\displaystyle {\displaystyle \sum_{j=0}^{min\left(a,b,M\right)}}\,\left[\begin{array}{c}
a\\
j
\end{array}\right]\left[\begin{array}{c}
b\\
M-j
\end{array}\right]}={\displaystyle \sum_{j=0}^{min\left(a,b\right)}}\left(\left[x^{j}\right]\left(x\right)_{a}\right)\cdot\left(\left[x^{M-j}\right]\left(x\right)_{b}\right)$$
$$=\left[x^{M}\right]\left(\left(x\right)_{a}\cdot\left(x\right)_{b}\right) \tag{2}\label{3}$$
Remark. The upper limit can be extended but additional terms would be zero.
The intent is to convert (1) to (2). Which can be done using:
http://functions.wolfram.com/06.10.17.0004.02
Substituting our symbols.
$$\left(x\right)_{b}=\frac{\Gamma\left(x+a\right)\Gamma\left(x+b\right)}{\Gamma\left(x\right)\Gamma\left(x+a+b\right)}\left(x+a\right)_{b}$$
Which is so obvious it's laughable.
Thus the answer is:
$$\left[x^{M}\right]\frac{\Gamma\left(x+a\right)\Gamma\left(x+b\right)}{\Gamma\left(x\right)\Gamma\left(x+a+b\right)}\left(x\right)_{\left(a+b\right)}=\left[x^{M}\right]\left(x\right)_{a}\cdot\left(x\right)_{b}$$
Maxima example code for all M:  
load ("stirling")$
gamma_expand: true$
gamma(x+5);
p1:pochhammer(x,6);
expand(p1);
p2:pochhammer(x+a,6);
conv(a,b):=pochhammer(x,a+b)*gamma(x+a)*gamma(x+b)/((gamma(x))*gamma(x+a+b));
ratsimp(conv(4,5));
ratsimp(pochhammer(x,4)*pochhammer(x,5));


