Estimate on sum of squares of multinomial coefficients I am interested in approximating the sum of the squares of the multinomial coefficients, i.e.
$a_\ell^p := \sum_{k_0+\ldots+k_p = \ell} (\frac{\ell!}{k_0! \ldots k_p!})^2$ 
or more general, 
$a_\ell^{\alpha_0,\ldots, \alpha_p} := \sum_{k_0+\ldots+k_p = \ell} (\prod_{i=0}^p \alpha_i^{k_i})^2(\frac{\ell!}{k_0! \ldots k_p!})^2$ 
Here $p$ is prime, and $\ell$ is an integer smaller then $p$, $k_i$ are non-negative integers. I would like to obtain some simple expression in terms of $\ell$ and $p$, a good approximation for large $p$, an upper bound will be good. 
I saw some results on recursive formula for such expressions, but not estimates. Should I just go with a Stirling formula in all the terms or is there something better done in this direction? 
Does someone knows what would be the Mathematica/Maple code for calculating such sums as functions from $p$ and $l$?
Thanks!
 A: Douglas already commented that the asymptotics for fixed $p$ and $l\to \infty$ shoudl follow from standard methods. One gets
$$a_{\ell}^p\approx (p+1)^{2\ell+\frac{p+1}{2}}(4\pi \ell)^{-\frac{p}{2}}.$$
See theorem 4 in "Counting Abelian squares", by Richmond and Shallit. Notice that these numbers appear also in combinatorics when considering abelian squares, or more generally abelian powers, on a fixed alphabet.
For the asymptotics that you're interested in, at least in the unweighted case, one can say
$$a _{\ell} ^p=\sum _{j=0} ^{\ell} \binom{p}{j}\sum _{a _1+ \cdots +a _j = \ell \atop a _i \geq 1} \binom{\ell}{a _1,a _2,\dots,a _j}^2$$
which makes it clear that $a _{\ell}^p$ is a polynomial in $p$ of fixed degree $\ell$. The coefficient of $\binom{p}{\ell}$ is $(\ell!)^2$, and the coefficient of $\binom{p}{\ell -1}$ is $\frac{\ell-1}{4}(\ell!)^2$, so you have
$$a _{\ell}^p =\ell!p^{\ell}-\ell!\frac{\ell(\ell-1)}{4}p^{\ell-1}+O(p^{\ell-2}).$$
A: Gjergji has given a very nice result, just want to explain Gjergji's answer for 
$$ a^p_l = l!p^l - l!\frac{l(l-1)}{4} p^{l-1} + O(p^{l-2})$$
To determine the coeffients of $p^l$, we see that this order is only included in $p \choose l$ because 
$$ {p\choose l} = \frac{(p-l+1)(p-l+2)\cdots (p)}{l!}$$
And in this case the term $\sum_{\substack{a_1 + \ldots + a_l = l \\ a_i \geq 1}} {l \choose {a_1, a_2, \ldots,a_l }}^2 = (l!)^2$ because we can only choose $a_i = 1, \forall i$. Therefore the coefficient of $p^l$ is $l!$. To determine the coefficient in $p^{l-1}$, we can see that it is contibuted by both $p \choose l$ and $ p \choose l-1$, for the term $p \choose l$, it contributes
$$ -(1 + 2 +\ldots + l-1)(l!) = -\frac{l(l-1)}{2}l!$$
For the term comtributed by $p \choose l-1$, notice $\sum_{\substack{a_1 + \ldots + a_{l-1} = l \\ a_i \geq 1}} {l \choose {a_1, a_2, \ldots,a_{l-1} }}^2 = (l-1)\frac{(l!)^2}{4}$, where the $(l-1)$ factor comes from the summation. Therefore the coefficient by $p \choose l-1$ is $(l-1)\frac{(l!)^2}{4} \frac{1}{(l-1)!} = l!\frac{l(l-1)}{4}$. So totally the coefficent is $-l!\frac{l(l-1)}{4}$, which is what Gjergji got in his answer. 
That's actually how the minus sign comes from.
