Drawing balls from an urn or counting certain posets

It would appear that these are multisets of multisets which may be enumerated using the Polya Enumeration Theorem (PET), same as what was posted at the following MSE link.

The answer is a special case of what was presented there and is given by

$$\left[\prod_{k=1}^q A_k^{N/q}\right] Z\left(S_{N/M}; Z\left(S_M; \sum_{k=1}^q A_{k}\right)\right).$$

In terms of combinatorial classes we have made use of the unlabeled class

$$\def\textsc#1{\dosc#1\csod} \def\dosc#1#2\csod{{\rm #1{\small #2}}} \textsc{MSET}_{=N/M} \left(\textsc{MSET}_{=M} \left(\sum_{k=1}^q \mathcal{A}_{k}\right)\right).$$

The linked post documents how to compute the cycle index of the symmetric group and how to substitute $Z(S_M)$ into $Z(S_{N/M}).$

The question is, why can we use multisets here? The answer is that there is a one-to-one mapping between those multisets and what OP calls building blocks. Every block obviously corresponds to exactly one multiset and every multiset to one block, namely by ordering its constituents. The same thing happens with multisets of blocks.

This was implemented in Maple and here are a few sample results that a potential contributor may compare to their work and / or use to verify that we have the correct understanding of the problem.

> seq(A(2,12,M), M in divisors(12));
                              1, 4, 5, 5, 4, 1

> seq(A(3,12,M), M in divisors(12));
                            1, 15, 25, 23, 10, 1

> seq(A(4,12,M), M in divisors(12));
                            1, 47, 93, 70, 22, 1

> seq(A(4,20,M), M in divisors(20));
                          1, 306, 2505, 1746, 73, 1

> seq(A(5,20,M), M in divisors(20));
                        1, 2021, 19834, 11131, 191, 1

The Maple code for the above goes as follows. The reader is invited to present their results in a language of their choice for a Rosetta stone effect.

with(combinat);
with(numtheory);

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;

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

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

    subsl := [];

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

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

        subsl := [op(subsl), v=subs(subs1, poly)];
    od;

    subs(subsl, ind);
end;

pet_cycleind_comp :=
proc(idxtrg, idx)
local varstrg, vars, vt, sbstrg, sbs, len;

    varstrg := indets(idxtrg);
    vars := indets(idx);

    sbstrg := [];

    for vt in varstrg do
        len := op(1, vt);

        sbs :=
        [seq(v = a[op(1, v)*len], v in vars)];

        sbstrg :=
        [op(sbstrg),
         a[len] = subs(sbs, idx)];
    od;

    expand(subs(sbstrg, idxtrg));
end;

A :=
proc(q, N, M)
option remember;
local cind_block, cind_word, cind_comp,
    vars, gf, vidx;

    if N mod q > 0 or N mod M > 0 then
        return FAIL;
    fi;

    cind_block := pet_cycleind_symm(M);
    cind_word := pet_cycleind_symm(N/M);

    cind_comp := pet_cycleind_comp(cind_word, cind_block);

    vars := add(A[p], p=1..q);
    gf := expand(pet_varinto_cind(vars, cind_comp));

    for vidx to q do
        gf := coeff(gf, A[vidx], N/q);
    od;

    gf;
end;

Addendum. With the above answer we compute the compound cycle index of the operator

$$\textsc{MSET}_{=N/M}(\textsc{MSET}_{=M}(\cdot))$$

and then substitute the variables into this cycle index. M. Scheuer in his answer suggests a different approach where we substitute the variables into the first cycle index, obtaining the multisets, and then substitute into $Z_{N/M}.$ Experimental data indicate that this approach is superior. The same effect can be achieved by not expanding the compound cycle index into its constitutents. This yields the following Maple code (duplicate code omitted).

pet_cycleind_comp :=
proc(idxtrg, idx)
local varstrg, vars, vt, sbstrg, sbs, len;

    varstrg := indets(idxtrg);
    vars := indets(idx);

    sbstrg := [];

    for vt in varstrg do
        len := op(1, vt);

        sbs :=
        [seq(v = a[op(1, v)*len], v in vars)];

        sbstrg :=
        [op(sbstrg),
         a[len] = subs(sbs, idx)];
    od;

    subs(sbstrg, idxtrg);
end;


A :=
proc(q, N, M)
option remember;
local cind_block, cind_word, cind_comp,
    vars, gf, vidx;

    if N mod q > 0 or N mod M > 0 then
        return FAIL;
    fi;

    cind_block := pet_cycleind_symm(M);
    cind_word := pet_cycleind_symm(N/M);

    cind_comp := pet_cycleind_comp(cind_word, cind_block);

    vars := add(A[p], p=1..q);
    gf := expand(pet_varinto_cind(vars, cind_comp));

    for vidx to q do
        gf := coeff(gf, A[vidx], N/q);
    od;

    gf;
end;

AX :=
proc(q, N, M)
option remember;
local cind_block, cind_word, cind_comp,
    vars, blocks, gf, vidx;

    if N mod q > 0 or N mod M > 0 then
        return FAIL;
    fi;

    cind_block := pet_cycleind_symm(M);
    vars := add(A[p], p=1..q);
    blocks := pet_varinto_cind(vars, cind_block);

    cind_word := pet_cycleind_symm(N/M);
    gf := expand(pet_varinto_cind(blocks, cind_word));

    for vidx to q do
        gf := coeff(gf, A[vidx], N/q);
    od;

    gf;
end;

This answer is a supplement to the great answer of @MarkoRiedel and a kind of reflection to better see the mechanisms working. The original problem was to determine $A_4(100,4)$ which can be written according to Markos answer as \begin{align*} A_4(100,4)=[a^{25}b^{25}c^{25}d^{25}]Z\left(S_{25}(Z\left(S_4;a+b+c+d\right)\right) \end{align*} It is a hopeless task to tackle this manually. But we can see all important aspects by taking a smaller parameter $N=8$ instead of $N=100$. Leaving the other parameter $M=4,q=4$ unchanged we have $N/M=N/q=2$ and we calculate \begin{align*} \color{blue}{A_2(8,4)=[a^{2}b^{2}c^{2}d^{2}]Z\left(S_{2}(Z\left(S_4;a+b+c+d\right)\right)}\tag{1} \end{align*}

Calculation of $Z(S_4)$:

We start calculating $Z(S_4)$ by using the recurrence relation: \begin{align*} Z(S_n)=\frac{1}{n}\sum_{j=1}^n a_j Z(S_{n-j})\qquad\textrm{where}\qquad Z(S_0)=1 \end{align*} We obtain \begin{align*} Z(S_0)&=1\\ Z(S_1)&=\frac{1}{1}\sum_{j=1}^1a_jZ(S_{1-j})=a_1Z(S_0)\\ &=a_1\\ Z(S_2)&=\frac{1}{2}\sum_{j=1}^2 a_jZ(S_{2-j})=\frac{1}{2}\left(a_1Z(S_1)+a_2Z(S_0)\right)\\ &=\frac{1}{2}\left(a_1^2+a_2\right)\tag{2}\\ Z(S_3)&=\frac{1}{3}\sum_{j=1}^3 a_jZ(S_{3-j})=\frac{1}{3}\left(a_1\frac{1}{2}\left(a_1^2+a^2\right)+a_2a_1+a_3\right)\\ &=\frac{1}{6}\left(a_1^3+3a_1a_2+2a_3\right)\\ \color{blue}{Z(S_4)}&=\frac{1}{4}\sum_{j=1}^4a_jZ(S_{4-j})\\ &=\frac{1}{4}\left(a_1\frac{1}{6}\left(a_1^3+3a_1a_2+2a_3\right)+a_2\frac{1}{2}\left(a_1^2+a_2\right)+a_3a_1+a_4\right)\\ &\,\,\color{blue}{=\frac{1}{24}\left(a_1^4+6a_1^2a_2+8a_1a_3+3a_2^2+6a_4\right)}\tag{3} \end{align*}

An ordinary generating function of the cycle index $Z(S(n))$ is \begin{align*} \sum_{n=0}^\infty Z(S_n)) z^n=\exp\left(a_1z+\frac{a_2}{2}z^2+\frac{a_3}{3}z^3+\cdots\right) \end{align*} We can use this function to make a plausibility check of (3) via Wolfram Alpha by typing

SeriesCoefficient[Exp[a_1*z+a_2*z^2/2+a_3*z^3/3+a_4*z^4/4],{z,0,4}]

which gives the expected result.

Calculation of $Z(S_4;a+b+c+d)$

We substitute $a+b+c+d$ in (2)

\begin{align*} Z(S_4)=\frac{1}{24}\left(a_1^4+6a_1^2a_2+8a_1a_3+3a_2^2+6a_4\right) \end{align*}

by replacing $a_j$ with $a^j+b^j+c^j+d^j$ and obtain \begin{align*} \color{blue}{Z(S_4;}&\color{blue}{a+b+c+d)}\\ &=\frac{1}{24}\left((a+b+c+d)^4+6(a+b+c+d)^2(a^2+b^2+c^2+d^2)\right.\\ &\qquad\qquad \left.+8(a+b+c+d)(a^3+b^3+c^3+d^3)\right.\\ &\qquad\qquad \left.+3(a^2+b^2+c^2+d^2)^2+6(a^4+b^4+c^4+d^4)\right)\\ &=\cdots\\ &\,\,\color{blue}{=a^4+a^3b+a^3c+a^3d+a^2b^2+a^2bc+a^2bd+a^2c^2+a^2cd+a^2d^2}\\ &\qquad\color{blue}{+ab^3+ab^2c+ab^2d+abc^2+abcd+abd^2+ac^3+ac^2d+acd^2+ad^3}\\ &\qquad\color{blue}{+b^4+b^3c+b^3d+b^2c^2+b^2cd+b^2d^2+bc^3+bc^2d+bcd^2+bd^3}\\ &\qquad\color{blue}{+c^4+c^3d+c^2d^2+cd^3+d^4}\tag{4} \end{align*}

Note that (4) has $35$ summands each with coefficient $1$, so that $Z(S_4;a+b+c+d)|_{a=b=c=d=1}=35$. This corresponds to (1) in the original post which is, since $M=4$ and $q=4$ \begin{align*} \sum_{1\leq x_1\leq x_2\leq x_3 \leq x_4\leq 4}1=\binom{4+4-1}{4}=\binom{7}{3}=35 \end{align*}

Calculation of $[a^2b^2c^2d^2]Z(S_2;Z(S_4;a+b+c+d))$:

We substitute (4) in $Z(S_2)=\frac{1}{2}\left(a_1^2+a_2\right)$ according to (2) and we obtain \begin{align*} \color{blue}{[a^2}&\color{blue}{b^2c^2d^2]Z(S_2;Z(S_4;a+b+c+d))}\\ &=[a^2b^2c^2d^2]\frac{1}{2}\left(a^4+a^3b+a^3c+a^3d+a^2b^2+a^2bc+a^2bd\right.\\ &\qquad\qquad\qquad\quad\left.+a^2c^2+a^2cd+a^2d^2+ab^3+ab^2c+ab^2d\right.\\ &\qquad\qquad\qquad\quad\left.+abc^2+abcd+abd^2+ac^3+ac^2d+acd^2\right.\\ &\qquad\qquad\qquad\quad\left.+ad^3+b^4+b^3c+b^3d+b^2c^2+b^2cd +b^2d^2 \right.\\ &\qquad\qquad\qquad\quad\left.+bc^3+bc^2d+bcd^2+bd^3+c^4+c^3d+c^2d^2\right.\\ &\qquad\qquad\qquad\quad\left.+cd^3+d^4\right)^2\\ &\quad+[a^2b^2c^2d^2]\frac{1}{2}\left(a^8+a^6b^2+a^6c^2+a^6d^2+a^4b^4+a^4b^2c^2+a^4b^2d^2\right.\\ &\quad\qquad\qquad\qquad\quad\left.+a^4c^4+a^4c^2d^2+a^4d^4 +a^2b^6+a^2b^4c^2+a^2b^4d^2\right.\\ &\quad\qquad\qquad\qquad\quad\left.+a^2b^2c^4+(abcd)^2+a^2b^2d^4+a^2c^6 +a^2c^4d^2+a^2c^2d^4 \right.\\ &\quad\qquad\qquad\qquad\quad\left.+a^2d^6+b^8+b^6c^2+b^6d^2+b^4c^4+b^4c^2d^2+b^4d^4\right.\\ &\quad\qquad\qquad\qquad\quad\left.+b^2c^6+b^2c^4d^2+b^2c^2d^4+b^2d^6 +c^8+c^6d^2+c^4d^4\right.\\ &\quad\qquad\qquad\qquad\quad\left.+c^2d^6+d^8\right)\\ &=[a^2b^2c^2d^2]\frac{1}{2}\left(a^2b^2+a^2bc+a^2bd+a^2c^2+a^2cd+a^2d^2\right.\\ &\qquad\qquad\qquad\quad\left.+ab^2c+ab^2d+abc^2+abcd+abd^2+ac^2d+acd^2\right.\\ &\qquad\qquad\qquad\quad\left.+b^2c^2+b^2cd+b^2d^2+bc^2d+bcd^2+c^2d^2\right)^2\\ &\quad+[a^2b^2c^2d^2]\frac{1}{2}(abcd)^2\tag{5}\\ &=[a^2b^2c^2d^2]\frac{1}{2}\left(2\left(a^2b^2\right)\left(c^2d^2\right) +2\left(a^2bc\right)\left(bcd^2\right) +2\left(a^2bd\right)\left(bc^2d\right)\right.\\ &\qquad\qquad\qquad\quad\left.+2\left(a^2c^2\right)\left(b^2d^2\right) +2\left(a^2cd\right)\left(b^2cd\right) +2\left(a^2d^2\right)\left(b^2c^2\right)\right.\\ &\qquad\qquad\qquad\quad\left.+2\left(ab^2c\right)\left(acd^2\right) +2\left(ab^2d\right)\left(ac^2d\right) +2\left(abc^2\right)\left(abd^2\right)\right.\\ &\qquad\qquad\qquad\quad\left.+(abcd)^2\right)\\ &\quad+[a^2b^2c^2d^2]\frac{1}{2}(abcd)^2\tag{6}\\ &\,\,\color{blue}{=10}\tag{7} \end{align*}

Comment:

• In (5) we keep only terms from the line above which have linear or square factors, since other terms do not contribute to $[a^2b^2c^2d^2]$.

• In (6) we multiply out and indicate the factors by keeping them in parenthesis.

We finally got the result (7): \begin{align*} \color{blue}{A_4(8,4)}&=[a^2b^2c^2d^2]Z(S_2;Z(S_4;a+b+c+d))\\ &\,\,\color{blue}{=10} \end{align*}

We can verify the result by listing the $35$ valid configurations according to $Z(S_4;a+b+c+d)$ $$ \begin{array}{cccc} 1111&1222&2222&3333\\ 1112&\color{blue}{1223}&2223&3334\\ 1113&\color{blue}{1224}&2224&\color{blue}{3344}\\ 1114&\color{blue}{1233}&\color{blue}{2233}&3444\\ \color{blue}{1122}&\color{blue}{1234}&\color{blue}{2234}&4444\\ \color{blue}{1123}&\color{blue}{1244}&\color{blue}{2244}&\\ \color{blue}{1124}&1333&2333&\\ \color{blue}{1133}&\color{blue}{1334}&\color{blue}{2334}&\\ \color{blue}{1134}&\color{blue}{1344}&\color{blue}{2344}&\\ \color{blue}{1144}&1444&2444&\\ \end{array} $$

The valid strings which have no more than two occurrences of each of the characters are marked $\mathrm{\color{blue}{blue}}$. Corresponding with $[a^2b^2c^2d^2]Z(S_2(Z(S_4;a+b+c+d))$ we concatenate the valid $\mathrm{\color{blue}{blue}}$ strings and find the $10$ resulting strings \begin{align*} 1122.3344\qquad1144.2233\\ 1123.2344\qquad1223.1344\\ 1124.2334\qquad1224.1334\\ 1133.2244\qquad1233.1244\\ 1134.2234\qquad1234.1234\\ \end{align*}