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!