Partitioning a multiset into multisets of fixed sizes Say we have a multiset $S(\mathbf{d}$) where $\mathbf{d}$ is a list of $l$ numbers and the multiplicity of the $i$th element of $S$ is $d_i$. The cardinality $N$ of $S$ is $\sum d_i$.
We want to partition $S$ into $m$ multisets of size $k_i$ respectively, so that $\sum k_i = \sum d_i = N$. How many ways can we do this?
In my mind this is a generalization of the multinomial coefficient $\binom{n}{k_1,k_2,\ldots,k_m}$ representing the number of ways to partition a set of $n=\sum k_i$ objects into $m$ bins of sizes $k_i$, to a sort of number like $\binom{\mathbf{d}}{k_1,k_2,\ldots,k_m}$ or $\binom{\mathbf{d}}{\mathbf{k}}$ representing the number of ways to partition a multiset of $n=\sum k_i = \sum d_i$ into $m$ bins of sizes $k_i$.
There are a few special cases that are simpler to calculate:


*

*If $m=1$, then clearly $k_1 = N$ and you're choosing the whole multiset. So $\binom{\mathbf{d}}{(N)} = 1$

*If $m=2$, then you only have to handle choosing $k_1$ or $k_2$ elements from a multiset, because the rest will be the other set. So, as mentioned below, you can use a generating function and $\binom{\mathbf{d}}{(k_1,k_2)}$ is equal to the coefficient of $x^{k_1}$ or $x^{k_2}$ in $\prod\limits_{i=1}^l 1 + x^2 + \cdots + x^{d_i} = \prod\limits_{i=1}^l \frac{1-x^{d_i - 1}}{1 - x}$. But then you also need to account for the fact that order doesn't matter, which I'm not sure how to do. Like in the first example below, you would find that there are $3$ ways to choose $2$ elements, but there are only $2$ ways to split the multiset because you have to choose 2 of them that are compatible.


Examples
Let's say that $\mathbf{d} = (2, 2)$, so $S(\mathbf{d})$ might be $\{a, a, b, b\}$. Let $k_1 = k_2 = 2$, so we need to find all the ways of splitting $S$ into two sub-multisets of size $2$. There are exactly $2$ ways of doing this: $\{\{a,a\},\{b,b\}\}$ and $\{\{a,b\},\{a,b\}\}$, so $\binom{(2,2)}{(2,2)} = 2$.
Another example: $\mathbf{d} = (2,2)$, so $S(\mathbf{d})$ could be $\{a,a,b,b\}$. Let $k_1 = 1$, $k_2 = 1$, and $k_3 = 2$. There are $3$ ways of doing this: $\{\{a\},\{a\},\{b,b\}\}$, $\{\{b\},\{b\},\{a,a\}\}$, and $\{\{a\},\{b\},\{a,b\}\}$. So $\binom{(2,2)}{(1,1,2)}=3$.
My Attempts
I've tried to figure this out two ways. The first was to find a recurrence relation and some base cases, kind of how Stirling numbers of the second kind can be computed using the identity $S(n,k) = kS(n-1,k) + S(n-1,k-1)$. I tried to think about what happens if you already have a partition and want to add an element to the original multiset, but then you have to decide which bin that element will go into or whether or not to add a new bin.
I also tried to derive it in the way that multinomial coefficients are derived, by counting the number of ways to fill the first bin, and then the second, and so on. The number of ways to choose $k_1$ elements from the multiset to put in the first bin can be computed by finding the coefficient of $x^{k_1}$ in $\prod\limits_{i=1}^l 1+x+x^2+\cdots+x^{d_i}$, which isn't explicit but it's a start. But then, depending on which elements you chose, you don't know how to adjust your multiset to reflect the remaining elements.
 A: I'm posting an implementation of Marko Riedel's algorithm in Sage because Sage is open source and widely available. It works on all the examples he posted, but for larger examples like $\binom{49, 49, 1, 1}{10, 10, 10, 10, 10, 10, 10, 10, 10, 10}$ it's hanging.
#!/usr/bin/env sage

import sys
from sage.all import *

Sym = SymmetricFunctions(QQ)
p = Sym.powersum()

def sub_cycle_index(Zout, Zin):
    """Substitutes Zin into Zout to produce Zout evaluated at Zin.

    This is accomplished by replacing p[i] in Zout with Zin, but with
    every p[j] in Zin replaced with p[i*j].
    """

    return p.sum(c * p.prod(Zin.frobenius(i) for i in mu) for mu, c in Zout)

def multiset_cycle_index(ms):
    """The cycle index of the given multiset, given by the formula

    .. math:: \prod_{\{k\}}\left( Z(S_{\sigma_k}; Z(S_k))\right)

    where :math:`\{k\}` are the elements of the multiset and
    :math:`\sigma_k` is the multiplicity of the element :math:`k`.
    """

    Z = lambda i: SymmetricGroup(i).cycle_index()
    return p.prod(sub_cycle_index(Z(sk), Z(k)) for k, sk in ms.items())

def list_to_multiset(els):
    """Converts a list of elements representing a multiset to a dictionary
    where the keys are the elements of the multiset and the values are
    the multiplicities.
    """

    ms = {}
    for x in els:
        ms[x] = ms.get(x,0) + 1
    return ms

def mset_choose(s, d):
    """Compute the "multiset coefficient" :math:`\binom{s}{d}`."""

    A = PolynomialRing(QQ, len(s), 'A').gens()
    mono = prod(a^i for a, i in zip(A, s))
    Z = multiset_cycle_index(list_to_multiset(d))
    return Z.expand(len(A), A).coefficient(mono)

if __name__ == '__main__':
    if len(sys.argv) != 3:
        print("Usage: %s 's_1, s_2, ..' 'd_1, s_2, ..'" % sys.argv[0])
        print("Outputs the number of ways the multiset s can be partitioned into multisets of sizes d_i.")
        sys.exit(1)

    s = map(int, sys.argv[1].split(' '))
    d = map(int, sys.argv[2].split(' '))

    if sum(s) != sum(d):
        print("The sum of the elements of s must equal the sum of the elements of d")
        sys.exit(1)

    print(mset_choose(s, d))

