If I'm not mistaken, you at least have the following bound:
$$
s(B) \ \; \leq\ \; d^n \,-\, \exp \left( -2 \left( \frac{\ln B}{\ln ( n! ) - d\ln ( (n/d)! )} \right)^2 \right) \,d^{n} .
$$
You may want to do a little fiddling around with Stirling's approximation to reduce the denominator in the fraction a bit further, but I didn't want to weaken the bound more than necessary.
I found this by upper-bounding the probability that a multinomial distribution with uniform bin probabilities yields an "atypical" observation, using Hoeffding's inequality. "Atypical" is here defined in the information-theoretic sense: An atypical observation is one whose logarithmic probability falls short of the expected logarithmic probability over the whole distribution (see Cover and Thomas' Elements of Information Theory, particularly chapter 11).
Some more details:
Let $p$ be point probabilities for a multinomial distribution with uniform bin probabilities:
$$
p(c) \ =\ \left( \frac{n!}{c_1!c_2!\cdots c_d!} \right) d^{-n}.
$$
Notice that
$$
\ln p(c) \,+\, n\ln d \,\ =\ \ln \frac{n!}{c_1!c_2!\cdots c_d!},
$$
and thus
$$
\ln p(c) \,+\, n\ln d \;\leq\; \ln B \;\quad \iff\quad \frac{n!}{c_1!c_2!\cdots c_d!} \;\leq\; B.
$$
Further, the entropy of the flat multinomial distribution is less than the entropy of $n$ samples from the flat categorical distribution with $d$ bins: The categorical includes order information as well as frequency information, while the multinomial only includes frequency information.
We thus have the following bound on the expected value of $-\ln p(c)$:
$$
E \left[\, -\ln p(c)\, \right] \ \leq\
n \cdot H\left( \frac{1}{d}, \frac{1}{d}, \ldots, \frac{1}{d} \right)
\ =\ n\ln d,
$$
or in other words, $-n\ln d \leq E[\ln p(c)]$.
Further, the minimum and maximum values of the random variable $\ln p(c)$ are
$$
a
\; =\;
\min_c \ln p(c)
\; =\;
\ln \frac{n!}{n!\, 0!\, 0!\, \cdots\, 0!} d^{-n}
\; =\;
-n\ln d;
$$
$$
b
\; =\;
\max_c \ln p(c)
\; =\;
\ln \frac{n!}{(n/d)!\, (n/d)!\, \cdots\, (n/d)!} d^{-n}
\; =\;
\ln ( n! ) - d\ln ( (n/d)! ) - n\ln d.
$$
The squared distance between these two extremes is consequently
$$
(a - b)^2 \ =\ \left( \ln ( n! ) - d\ln ( (n/d)! ) \right)^2.
$$
We can now make the following comparisons:
\begin{eqnarray}
\Pr\left\{ \frac{n!}{c_1!c_2!\cdots c_d!} < B \right\}
& = &
\Pr\left\{ \ln p(c) \,+\, n\ln d < \ln B \right\}
\\
& \leq &
\Pr\left\{ \ln p(c) \,-\, E[\ln p(c)] < \ln B \right\}
\\
& = &
1 \,-\, \Pr\left\{ \ln p(c) \,-\, E[\ln p(c)] \geq \ln B \right\}
\\
& \leq &
1 \,-\, \exp \left( -2 \left( \frac{\ln B}{\ln ( n! ) - d\ln ( (n/d)! )} \right)^2 \right).
\end{eqnarray}
The first inequality follows from the fact that a proposition always has a lower probability than its logical consequences; the second is the application of Heoffding's inequality.
We thus have
$$
\Pr\left\{ \frac{n!}{c_1!c_2!\cdots c_d!} < B \right\}
\ \leq \
1\;-\;\exp \left( -2 \left( \frac{\ln B}{\ln ( n! ) - d\ln ( (n/d)! )} \right)^2 \right).
$$
By multiplying by $d^n$, we obtain the inequality stated above, since the probability of a sample from a flat multinomial distribution is equal to the corresponding multinomial coefficient divided by $d^n$.
First, what the Stirling bound or Stanica's result give is already a $(1+O(n^{-1}))$ approximation of $\binom nk$, hence the only problem can be with the sum. I don't know how to do that with such precision, but it's easy to compute it up to a constant factor by approximating with a geometric series:
$$\sum_{i\le k}\binom ni=\begin{cases}\Theta(2^n)&k\ge n/2-\sqrt n,\\\\\Theta\left(\left(1-\frac{2k}n\right)^{-1}\binom nk\right)&k\le n/2-\sqrt n.\end{cases}$$
More generally,
$$\sum_{i\le k}\binom ni(1-\epsilon)^{n-i}\epsilon^i=\begin{cases}\Theta(1)&k\ge\epsilon n-s,\\\\\Theta\left(\frac{\epsilon(n-k)}{\epsilon n-k}\binom nk(1-\epsilon)^{n-k}\epsilon^k\right)&k\le\epsilon n-s,\end{cases}$$
where $s=\sqrt{n\epsilon(1-\epsilon)}$. Cf. the appendix to my paper http://math.cas.cz/~jerabek/papers/wphp.pdf .
Best Answer
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}).$$