Distributing groups of objects into boxes

The problem can viewed as enumerating $r\times k$ matrices with non-negative integer entries that have (i) row sums $a_1,\dots,a_r$, and (ii) pairwise distinct columns. The number of such matrices then needs to be divided by $k!$ since the order of columns (boxes) does not matter.

For given ranges of parameters the best shot seems to be using the inclusion-exclusion principle to account for condition (ii). It amounts to summing over the partitions of $k$ in the following formula: $$\frac{(-1)^k}{k!}\sum_{p\ \vdash\ k} (-1)^{|p|} \pi(p) \prod_{j=1}^r f(a_j,p),$$ where $\pi(p)$ is the number of permutations of cycle type $p$, and $f(a_j,p)$ for $p=(p_1,\dots,p_m)$ is the number of non-negative integer solutions $(n_1,\dots,n_m)$ to the equation: $$p_1n_1 + p_2n_2 + \dots + p_mn_m = a_j.$$

For each fixed partition $p=(p_1,\dots,p_m)$, the values $f(a_j,p)$ can be obtained at once for all $j\in\{1,\dots,r\}$ by the dynamic-programming algorithm based on the identity: \begin{split} f(a,p) &= [x^a]\ \frac1{(1-x)^{e_1}\cdots(1-x^k)^{e_k}} \\ &= \sum_{t=0}^{\lfloor a/k\rfloor} (-1)^t \binom{-e_k}{t} [x^{a-tk}] \frac1{(1-x)^{e_1}\cdots(1-x^{k-1})^{e_{k-1}}}, \end{split} where $e_i$ is the multiplicity of part $i$ in $p$, with running time $O(k\cdot \max_j a_j)$.

In the given example we have $k=4$, $r=2$ and $(a_1,a_2)=(4,2)$, and we need to sum up over the partitions of $k$, i.e., $$\{(1,1,1,1),(2,1,1),(3,1),(2,2),(4)\}.$$ Correspondingly, the summands are