Stirling numbers of the second kind on Multiset
There is no known formulation for a general multiset. However, a paper at JIS tackles the case where the element 1 occurs multiple times.
By way of enrichment here is a solution using Power Group Enumeration in polynomial form. This is not the most efficient algorithm and if we are just after the count we can in fact do better by enumerating all such set partitions recursively, which is also included here. The PGE algorithm has a certain chance of actually leading to a closed form expression rather than the output of some procedure, however.
The algorithm is the same as presented at the following MSE link I and at this link MSE link II (edge colorings of the cube). We require the cycle indices of the permutation groups acting on the slots and on the values being distributed into them. For a multiset
$$A_1^{q_1} A_2^{q_2} \times \cdots \times A_p^{q_p}$$
with $q_k$ the multiplicities the slots consist of $Q= \sum_r q_r$ distinct slots, into which we distribute the identifiers of the multisets where that slot is to reside. We may imagine $q_r$ slots of type $r$ forming a block, with the blocks being adjacent and containing only slots of the same type. Now with a block forming a multiset the group acting on these slots is the symmetric group $S_{q_r}$ with cycle index $Z(S_{q_r}).$ Therefore the cycle index of the slot permutation group (all blocks) is $$\prod_r Z(S_{q_r}).$$ On the other hand the multiset identifiers are permuted by the symmetric group $S_Q$ with cycle index $Z(S_Q).$ These data suffice to run the algorithm. Observe that when we cover (consult links for the details of this procedure) the cycles of a slot permutation with a cycle from a multiset identifier permutation we must record the identifiers used in the covering in a generating function. We then obtain the desired result by replacing the indeterminates in the terms of this generating function (once computed) with $v^m$ where $m$ is the number of indeterminates in the term. This yields the generating function by the number of multisets for all $k$ at once.
We present some examples. When all $q_r$ are one we are working with $p$ distinct items and should obtain Stirling numbers. (Stirling numbers exemplify the complexity of this problem as the PGE routine is faster here and takes less memory than the recursive enumeration routine. This is because the slot permutation group in this case contains only one element. The computation was done at this MSE link III.) Indeed we find
> add(v^k*stirling2(3,k), k=1..3);
3 2
v + 3 v + v
> MSETS([1,1,1]);
3 2
v + 3 v + v
> add(v^k*stirling2(7,k), k=1..7);
7 6 5 4 3 2
v + 21 v + 140 v + 350 v + 301 v + 63 v + v
> MSETS([seq(1, w=1..7)]);
7 6 5 4 3 2
v + 21 v + 140 v + 350 v + 301 v + 63 v + v
> add(v^k*stirling2(10,k), k=1..10);
10 9 8 7 6 5
v + 45 v + 750 v + 5880 v + 22827 v + 42525 v
4 3 2
+ 34105 v + 9330 v + 511 v + v
> MSETS([seq(1, w=1..10)]);
10 9 8 7 6 5
v + 45 v + 750 v + 5880 v + 22827 v + 42525 v
4 3 2
+ 34105 v + 9330 v + 511 v + v
On the other hand when $p=1$ we have $q_1$ like items and should get the values of the partition function $p_k(n).$ Here we find once again
> add(v^k*part(3,k), k=1..3);
3 2
v + v + v
> MSETS([3]);
3 2
v + v + v
> add(v^k*part(7,k), k=1..7);
7 6 5 4 3 2
v + v + 2 v + 3 v + 4 v + 3 v + v
> MSETS([7]);
7 6 5 4 3 2
v + v + 2 v + 3 v + 4 v + 3 v + v
> add(v^k*part(10,k), k=1..10);
10 9 8 7 6 5 4 3 2
v + v + 2 v + 3 v + 5 v + 7 v + 9 v + 8 v + 5 v + v
> MSETS([10]);
10 9 8 7 6 5 4 3 2
v + v + 2 v + 3 v + 5 v + 7 v + 9 v + 8 v + 5 v + v
Now of course it gets difficult and interesting when we have less structure in the size and number of the source multiset. Here are a few last examples:
> ENUM([2,3,4]);
9 8 7 6 5 4 3 2
v + 6 v + 26 v + 75 v + 151 v + 191 v + 126 v + 29 v
+ v
> MSETS([2,3,4]);
9 8 7 6 5 4 3 2
v + 6 v + 26 v + 75 v + 151 v + 191 v + 126 v + 29 v
+ v
> ENUM([3,4,5]);
12 11 10 9 8 7 6
v + 6 v + 30 v + 111 v + 328 v + 757 v + 1331 v
5 4 3 2
+ 1652 v + 1264 v + 474 v + 59 v + v
> MSETS([3,4,5]);
12 11 10 9 8 7 6
v + 6 v + 30 v + 111 v + 328 v + 757 v + 1331 v
5 4 3 2
+ 1652 v + 1264 v + 474 v + 59 v + v
> ENUM([2,2,3,3]);
10 9 8 7 6 5 4
v + 10 v + 63 v + 257 v + 704 v + 1223 v + 1204 v
3 2
+ 536 v + 71 v + v
> MSETS([2,2,3,3]);
10 9 8 7 6 5 4
v + 10 v + 63 v + 257 v + 704 v + 1223 v + 1204 v
3 2
+ 536 v + 71 v + v
The source code is quite compact, with the answer being about half and the rest for testing and verification.
with(combinat);
pet_cycleind_symm :=
proc(n)
option remember;
if n=0 then return 1; fi;
expand(1/n*add(a[l]*pet_cycleind_symm(n-l), l=1..n));
end;
MSETS :=
proc(mset)
option remember;
local idx_slots, idx_sets, res, term_a, term_b,
v_a, v_b, inst_a, inst_b, len_a, len_b, p, q;
idx_slots :=
expand(mul(pet_cycleind_symm(item), item in mset));
if not(type(idx_slots, `+`)) then
idx_slots := [idx_slots];
fi;
idx_sets :=
pet_cycleind_symm(add(item, item in mset));
if not(type(idx_sets, `+`)) then
idx_sets := [idx_sets];
fi;
res := 0;
for term_a in idx_slots do
for term_b in idx_sets do
p := 1;
for v_a in indets(term_a) do
len_a := op(1, v_a);
inst_a := degree(term_a, v_a);
q := 0;
for v_b in indets(term_b) do
len_b := op(1, v_b);
inst_b := degree(term_b, v_b);
if len_a mod len_b = 0 then
q := q +
len_b*
add(mul(B[len_b,
len_b*cyc_num+cyc_pos],
cyc_pos=1..len_b),
cyc_num=0..inst_b-1);
fi;
od;
p := p*q^inst_a;
od;
res := res +
lcoeff(term_a)*lcoeff(term_b)*p;
od;
od;
map(term -> lcoeff(term)*v^nops(indets(term)),
expand(res));
end;
part :=
proc(n, k)
option remember;
local res;
if n=0 or k=0 then
return `if`(n=0 and k=0, 1, 0);
fi;
if k=1 then return 1 fi;
res := 0;
if n >= 1 then
res := res + part(n-1, k-1);
fi;
if n >= k then
res := res + part(n-k, k);
fi;
res;
end;
ENUM :=
proc(mset)
option remember;
local total, flat, recurse, msets, res;
total := add(item, item in mset);
flat :=
[seq(seq(p, q=1..mset[p]), p=1..nops(mset))];
msets := table(); res := 0;
recurse :=
proc(idx, len, part)
local pos, cur, sorted;
if idx > total then
sorted :=
sort([seq(part[pos], pos=1..len)]);
if not(type(msets[sorted], `integer`)) then
res := res + v^len;
msets[sorted] := 1;
fi;
return;
fi;
cur := A[flat[idx]];
for pos to len do
part[pos] := part[pos]*cur;
recurse(idx+1, len, part);
part[pos] := part[pos]/cur;
od;
part[len+1] := part[len+1]*cur;
recurse(idx+1, len+1, part);
part[len+1] := part[len+1]/cur;
end;
recurse(1, 0, Array([seq(1, q=1..total)]));
res;
end;
Here are two links to get you started:
Eulerian numbers of the second kind may be helpful (for counting ascents, descents, etc., though i think)
Additionally some more useful information may be found in Stanley's book