Proving a binomial identity involving the sum of a product of four binomial coefficients While working on a problem, I needed to prove the following identity:

For all $i,j,m,n \geqslant 0$,
  \begin{align*}
\sum_{k+\ell=m} k! \ell! \binom{n+k-1}{k} \binom{j-i}{k} \binom{n-j+m-1}{\ell} \binom{i}{\ell} = m! \binom{n-i+m-1}{m}\binom{j}{m}.
\end{align*}

Computational evidence indicates that this is true, but I haven't been able to prove it. It looks a little like Vandermonde's identity, but I haven't been able to massage it into a form where I can just apply that. Is there a clever way to prove this identity?
 A: Let 
\begin{align}
F(m,k)
&=k! (m-k)! \binom{n+k-1}{k} \binom{j-i}{k} \binom{n-j+m-1}{m-k} \binom{i}{m-k}
\\\\
&=\frac{(n+k-1)!}{(n-1)!}\binom{j-i}{k}\frac{(n-j+m-1)!}{(n-j+k-1)!}\binom{i}{m-k}
\end{align}
denote the summand. Let 
$$
S_m=\sum_{k=0}^m F(m,k)
$$ be the sum to be computed.
Define $G(m,k)$ as follows. The motivation for this definition is as of yet unclear$^{**}$:
\begin{align}
G(m,k)
&=F(m,k)\cdot \frac{k\,\left(m-k-i\right)\,\left(n+k-j-1\right)}
{m-k+1}
\\\\
&=\frac{(n+k-1)!}{(n-1)!}\binom{j-i}{k}\frac{(n-j+m-1)!}{(n-j+k-1)!}\color{red}{\binom{i+1}{m-k+1}}\frac{k\,\left(m-k-i\right)\,\left(n+k-j-1\right)}
{\color{red}{i+1}}
\end{align}
Note that the first expression for $G(m,k)$ is undefined when $k=m+1$, as it contains $\frac00$. However, it can be rearranged to the second expression, where $G(m,k)$ is unambiguously defined for all $k$. Note that $G(m,k)=0$ for all $k$ outside $[1,m+1]$. 
You can verify via tedious algebraic manipulations$^*$ that
$$
(m+1)F(m+1,k)-(j-m)(n+m-i)F(m,k)=G(m,k+1)-G(m,k)\label1\tag1
$$
holds for all nonnegative integer inputs. 
Now, sum both sides of the above equation from $k=0$ to $m+1$. The right had side is a telescoping sum. We attain
$$
(m+1)S_{m+1}-(j-m)(n+m-i)S_m=G(m,m+2)-G(m,0)=0-0,
$$
or
$$
S_{m+1}=\frac{(j-m)(n+m-i)}{m+1}S_m
$$
This easily lets you prove that $S_m=m!\binom{n-i+m-1}{m}\binom{j}m$ by induction.

$^*$ The calculations really are tedious, but certainly feasible for even a high schooler to do. Start by dividing $\eqref1$ though by $F(m,k)$, then replace all of the binomial coefficients with factorials. Every factorial term in the numerator will have a counterpart which cancels nicely with one in the denominator. For example, in $F(m+1,k)/F(m,k)$, there will appear $$\binom{i}{m+1-k}\Big/\binom{i}{m-k}=\frac{i!}{(m+1-k)!(i-m-1+k)!}\Big/\frac{i!}{(m-k)!(i-m+k)!}=\frac{i-m+k}{m+1-k}$$
After making all these cancellations, you will be left with a sum of fractions of polynomials in $i,j,k,m,n$. After clearing denominators, this is a polynomial equation which can be verified by collecting all like coefficients.
But there is not reason not to just have Mathematica verify the equation for you automatically.
F[m_,k_] := k! * (m-k)! * 
            Binomial[n+k-1,k] * Binomial[j-i,k] * 
            Binomial[n-j+m-1,m-k] * Binomial[i,m-k];

G[m_,k_] := F[m,k] * k * (m - k - i) * (n + k - j - 1) / (m - k + 1);

Print[
  FullSimplify[
    FunctionExpand[
    (m + 1)F[m + 1,k] - (j - m)(n + m - i)F[m,k] == G[m,k + 1] - G[m,k]
    ]
  ]
]

Try it online!

$^{**}$ The quotient $ \frac{k\,\left(m-k-i\right)\,\left(n+k-j-1\right)}{m-k+1}$ in the definition of $G$, and the coefficients $m+1$ and $i(j-m)(n+m-i)$ on the left hand side of $\eqref1$, are attained using Zeilberger's algorithm. Specifically, I used the programming language Maxima, this code snippet:
load (zeilberger)$

Zeilberger  
(
    k! * (m-k)! * 
    binomial(n+k-1,k) * binomial(j-i,k) * binomial(n-j+m-1,m-k) * binomial(i,m-k)
    ,k,m
);

Try it online! (ignore the small differences in the linked code).
This algorithm, and others is discussed in the book A = B by Marko Petkovsek, Herbert Wilf and Doron Zeilberger, available for free online. Basically, the book shows must summations similar to yours can be solved completely automatically, and computer verifiable proofs of these equation can also be produced automatically. 
A: \begin{align*}
\sum_{k+\ell=m} k! \ell! \binom{n+k-1}{k} \binom{j-i}{k} \binom{n-j+m-1}{\ell} \binom{i}{\ell} = m! \binom{n-i+m-1}{m}\binom{j}{m}.
\end{align*}
Okay, first thing I'd do, is expand the sum termwise a bit to cancel things to be:
\begin{align*}
\sum_{k+\ell=m} \frac{(n+k-1)!i!}{(n-1)!(i-\ell)!} \binom{j-i}{k} \binom{n-j+m-1}{\ell} = m! \binom{n-i+m-1}{m}\binom{j}{m}.
\end{align*}
Scanning a bit more, we see the possibility to get rid of $m!$ on the right hand side:
\begin{align*}
\sum_{k+\ell=m} \frac{(n+k-1)!i!}{(n-1)!(i-\ell)!} \binom{j-i}{k} \binom{n-j+m-1}{\ell} = \frac{(n-i+m-1)!}{(n-i-1)!}\binom{j}{m}.
\end{align*}
Using the fact (finally...) that $k+\ell=m$, noting that means $m!$ has $\ell!$ and $k!$ as factors, remaining factor p,we get:
\begin{align*}
\sum_{k+\ell=m} \frac{p(n+k-1)!(i)!(j-i)!(n-j+m-1)!}{(n-1)!(i-\ell)!(j-i-k)!(n-j+k-1)!} = \frac{p(n-i+m-1)!(j)!}{(n-i-1)!(j-m)!}.
\end{align*}
Looking a bit closer, we have $$\ell\leq i\leq j\leq n, m\leq n,k\leq j, k\leq n$$ As well as simplifying because factorials in denominators have a factor separating them, we can go down to:
\begin{align*}
\sum_{k+\ell=m} qsrt(n+k-1)!= uvp^2
\end{align*}
Where p,q,r,s,t,u, and v are all quotients of factorials. So, you have the statement that, a factorial and a square of a quotient of factorials, can be related by another set of 6 not necessarily distinct quotients of factorials, 
