Number of ways to distribute $n$ identical balls into $k$ identical boxes such that no box contains more than $m$ balls?

Question

I wonder how to count the number of ways (algorithmically is fine) to distribute $n$ identical balls into $k$ identical boxes such that no box contains more than $m$ balls?

I've run into answers in which either the balls are non-identical or the boxes are non-identical, but not this version.

Answer 1

This type of questions can be solved using Gaussian binomial coefficient.When we look at the its combinatorial usage , we can see that "Let $B(n,m,r)$ be the number of ways of throwing $r$ indistinguishable balls into $m$ indistinguishable bins, where each bin can contain up to $n$ balls. The Gaussian binomial coefficient can be used to characterize $B(n,m,r)$."

$$B(n,m,r)=[q^r]\binom{n+m}{m}_q$$ where $[q^r]P$ denotes the coefficient of $q^r$ in polynomial $P$

Answer 2

This can also be done using the Polya Enumeration Theorem. We get from first principles that the answer is given by where we use the cycle index of the symmetric group

$$Q_{n,k,m} = [A^n] Z(S_k; 1+A+A^2+\cdots+A^m).$$

We have the recurrence by Lovasz for the cycle index which says that

$$Z(S_k) = \frac{1}{k} \sum_{\ell=1}^k a_\ell Z(S_{k-\ell}) \quad\text{where}\quad Z(S_0) = 1.$$

Extracting coefficents we have the memoized recurrence

$$Q_{n,k,m} = [A^n] \frac{1}{k} \sum_{\ell=1}^k \sum_{p=0}^m A^{p\ell} Z(S_{k-\ell}; 1+\cdots+A^m) \\ = \frac{1}{k} \sum_{\ell=1}^k \sum_{p=0}^{\min(m,\lfloor n/\ell \rfloor)} [A^n] A^{p\ell} Z(S_{k-\ell}; 1+\cdots+A^m) \\ = \frac{1}{k} \sum_{\ell=1}^k \sum_{p=0}^{\min(m,\lfloor n/\ell \rfloor)} Q_{n-p\ell, k-\ell, m}$$

where the base cases are $Q_{n,0,m} = \delta_{0,n}$ and $Q_{n,1,m} = [[n\le m]].$ For example wth $k=5$ and $m=5$ we obtain the following sequence

$$1, 1, 2, 3, 5, 7, 9, 11, 14, 16, 18, 19, 20, 20, 19, 18, \\ 16, 14, 11, 9, 7, 5, 3, 2, 1, 1, 0, 0, 0, 0, 0, \ldots$$

which points us to OEIS A102422 where we see that we have the right data.

The code for this was as follows:

with(combinat);

pet_varinto_cind :=
proc(poly, ind)
local subs1, subs2, polyvars, indvars, v, pot, res, k;

    res := ind;

    polyvars := indets(poly);
    indvars := indets(ind);

    for v in indvars do
        pot := op(1, v);

        subs1 :=
        [seq(polyvars[k]=polyvars[k]^pot,
             k=1..nops(polyvars))];

        subs2 := [v=subs(subs1, poly)];

        res := subs(subs2, res);
    od;

    res;
end;

pet_cycleind_symm :=
proc(n)
local l;
option remember;

    if n=0 then return 1; fi;

    expand(1/n*add(a[l]*pet_cycleind_symm(n-l), l=1..n));
end;

X :=
proc(n,k,m)
local rep, gf, p;
option remember;
    
    rep := add(A^p, p=0..m);
    gf := pet_varinto_cind(rep, pet_cycleind_symm(k));

    coeff(expand(gf), A, n);
end;

Q :=
proc(n,k,m)
local l, p;
option remember;
    if k=0 then return `if`(n = 0, 1, 0) fi;
    if k=1 then return `if`(n <= m, 1, 0) fi;

    1/k*add(add(Q(n-p*l,k-l,m),
                p=0..min(m,floor(n/l))), l=1..k);
end;