Stirling numbers of the second kind on Multiset

Question

Stirling numbers of the second kind $S(n, k)$ count the number of ways to partition a set of $n$ elements into $k$ nonempty subsets. What if there were duplicate elements in the set? That is, the set is a multiset?

Answer 1

There is no known formulation for a general multiset. However, a paper at JIS tackles the case where the element 1 occurs multiple times.

Answer 2

Here are two links to get you started:

1. Eulerian numbers of the second kind may be helpful (for counting ascents, descents, etc., though i think)

2. Additionally some more useful information may be found in Stanley's book

Answer 3

You divide the Stirling number by the possible permutations of identical elements as this is just changing order within the stars and bars example.

Answer 4

For those of you who want a faster program to compute this, I have written a Matlab program for this. Since the code does not require symbolic math, it is a lot faster than Marko Reidel's Maple code, especially when the multiset is large. In addition, the algorithm does not require the enumeration of the set partitions but instead it is based on a recurrence relation given in Gupta (1961), "On the Partition of J-partite Numbers," Proceedings of the National Institute of Sciences of India Part A 27, 579-587, but I have added quite a few improvements in my code. Here is the code.

%
%   multipartr.m   Date: 12/19/2024
%   This Matlab function computes the number of
%   set partitions for (n1*1, n2*2, ..., ns*s) with 
%   exactly r parts for r=1,...,n1+...+ns
%   Input:
%   n: a vector that indicates the frequencies of
%      the multiset
%   Output:
%   y : a vector of the number of set partitions of 
%       the multiset that have exactly r parts for 
%       r=1,...,n1+...+ns
%
function y = multipartr(n)
   s = length(n);
   if s==1
      F = ones(n,'uint64');
      for j=2:n
          F(:,j) = F(:,j-1);
          F(j,j) = F(j,j)+1;
          for i=j+1:n
              F(i,j) = F(i,j)+F(i-j,j);
          end
      end
      y = [1 diff(F(n,:))]';
      return
   end
   n = sort(n,'descend');
   n1 = n+1;
   p = cumprod([1 n1(1:s-1)])';
   pn = prod(n1);
%
%  Create an index vector for \kappa's,
%  where 0<=\kappa<=n.  The index
%  maps each \kappa to an element
%  in the set of \kappa that is
%  descending.
%
   ind = zeros(pn,1);
   kappa = zeros(1,s);
   count = 1;
   sk = 0;
   for i=2:pn
       j = 1;
       while kappa(j)==n(j)
          kappa(j) = 0;
          sk = sk-n(j);
          j = j+1;
       end
       kappa(j) = kappa(j)+1;
       sk = sk+1;
       if all(diff(kappa)<=0)
          ind(i) = count;
          count = count+sk;
       else
          i1 = sort(kappa,'descend')*p+1;
          ind(i) = ind(i1);
       end
   end
   A = zeros(count,1,'uint64');
   A(1) = 1;
   kappa = zeros(1,s);
   sn = sum(n);
   sk = 0;   % sk = |\kappa|
%  loop through \kappa with 0<\kappa<=n
%  but only for those kappa's that are
%  descending.
   while sk<sn
      j = 1;
      while kappa(j)==n(j)
         j = j+1;
      end
      kappa(j) = kappa(j)+1;
      if j==1
         sk = sk+1;
      else
         kappa(1:j-1) = kappa(j);
         sk = sum(kappa);
      end
      ind1 = kappa*p+1;
      indc = ind(ind1);
      y1 = zeros(sk,1,'uint64');
      y1(1) = 1;
      nu = zeros(1,s);
      snu = 0;
      ind2 = 0;     % ind2 = nu*p
%  loop through \nu with 0<\nu<\kappa
      for ii=1:prod(kappa+1)-2
          j1 = 1;
          while nu(j1)==kappa(j1)
             nu(j1) = 0;
             snu = snu-kappa(j1);
             ind2 = ind2-kappa(j1)*p(j1);
             j1 = j1+1;
          end
          nu(j1) = nu(j1)+1;
          snu = snu+1;
          ind2 = ind2+p(j1);
          sk2 = sk;
          for r=1:min(floor(kappa./nu))
              sk2 = sk2-snu;
              if sk2>0
                 c1 = ind(ind1-r*ind2)-r;
                 for i=r+1:r+sk2
                     y1(i) = y1(i)+A(c1+i);
                 end
              else
                 y1(r) = y1(r)+1;
              end
          end
      end
      for i=1:sk
          A(indc+i) = y1(i)/i;
      end
   end
   y = A(end-sn+1:end);

If you do multipartr([2 2 3 3]), you get

      1
     71
    536
   1204
   1223
    704
    257
     63
     10
      1

Even for something as large as multipartr([12 12 13 13]), it takes 1.1 seconds to do on my Ryzen 5950x. If you are only interested in obtaining the total number of set partitions, i.e., sum(multipartr([12 12 13 13])), it can be done in 0.1 second. E-mail me if you are interested in having this code.

Answer 5

By definition, a set only has distinct elements. If elements are duplicated, then yes, the set in question would be a multiset.