Number of ways to distribute objects, some identical and others not, into identical groups

Using the notation from the following MSE link I we start with the source multiset

$$\def\textsc#1{\dosc#1\csod} \def\dosc#1#2\csod{{\rm #1{\small #2}}} \prod_{k=1}^l A_k^{\tau_k}$$

where we have $l$ different values and their multiplicities are the $\tau_k.$ We ask about the distinct partitions of this multiset into $N$ factors including one as a factor, where distinct refers to permutations of the $N$ factors by the symmetric group. The case where one is not admitted as a factor was discussed at the following MSE link II.

If we have a CAS like Maple, $N$ is reasonable and we seek fairly instant computation of these values then we may just use the cycle index $Z(S_N)$ of the symmetric group which implements the unlabeled operator $\textsc{MSET}_{=N}.$ This yields the formula

$$\left[\prod_{k=1}^l A_k^{\tau_k}\right] Z\left(S_N; \prod_{k=1}^l \frac{1}{1-A_k}\right).$$

Here we use the recurrence by Lovasz for the cycle index $Z(S_N)$, which is

$$Z(S_N) = \frac{1}{N} \sum_{l=1}^N a_l Z(S_{N-l}) \quad\text{where}\quad Z(S_0) = 1.$$

Maple can extract these coefficients by asking for the coefficient of the corresponding Taylor series. We get the following transcript:

> FACTORS(60^3, 3);
                                475

> FACTORS(60^4, 3);
                               1710

> FACTORS(120, 4); 
                                20

> FACTORS(512, 4);
                                18

> FACTORS(729, 5);
                                10

> FACTORS(2^4*3^3*5^2, 6);
                                573

> seq(FACTORS(n,4), n=1..60);
1, 1, 1, 2, 1, 2, 1, 3, 2, 2, 1, 4, 1, 2, 2, 5, 1, 4, 1, 4, 2, 2,

    1, 7, 2, 2, 3, 4, 1, 5, 1, 6, 2, 2, 2, 9, 1, 2, 2, 7, 1, 5, 1,

    4, 4, 2, 1, 11, 2, 4, 2, 4, 1, 7, 2, 7, 2, 2, 1, 11

The sequence is OEIS A218320 and looks to have the right values. The Maple code here is quite simple.

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, subs2, polyvars, indvars, v, pot, res;

    res := ind;

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

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

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

        subs2 := [v=subs(subs1, poly)];

        res := subs(subs2, res);
    od;

    res;
end;


MSETS :=
proc(src, N)
local msetgf, cind, gf, cf;

    msetgf := mul(1/(1-A[q]), q=1..nops(src));
    cind := pet_cycleind_symm(N);

    gf := pet_varinto_cind(msetgf, cind);

    for cf to nops(src) do
        gf := coeftayl(gf, A[cf] = 0, src[cf]);
    od;

    gf;
end;

FACTORS :=
proc(n, N)
local mults;

    mults := map(el -> el[2], op(2, ifactors(n)));
    MSETS(mults, N);
end;

Remark. An observation. While Burnside and Polya certainly represent an enrichment here we must also take care to include the basics, which in this case consists of a simple recurrence that is given at the OEIS entry and which computes the desired values almost instantly. With the variables renamed to indicate semantics we have an algorithm whose correctness follows by inspection and which is shown below.

FACTREC :=
proc(val, numel, maxfact)
option remember;
local divs;

    if numel = 1 then
        return `if`(val <= maxfact, 1, 0);
    fi;

    divs := select(d -> d <= maxfact, divisors(val));
    add(FACTREC(val/d, numel-1, d), d in divs);
end;

FACTORSX := (n, N) -> FACTREC(n, N, n);

This problem can be solved by an application of Burnside's lemma.

Let $X = \{(x,y,z) \in \mathbb N^3 : xyz = 60^3\}$ be the set of all ways to specify the cuboid where the order of the sides does matter. The group $G = S_3$ acts on the elements of $X$ by permuting the ordered triple $(x,y,z)$. We are looking for the number of orbits $|X/G|$.

To do this, we compute the number of elements of $X$ fixed by each element of $G$ and average them:

• $X^e$, the set of elements fixed by the identity element $e$, is just $X$. We have $|X| = |X^e| = \binom82 \binom52 \binom52$ by applying stars and bars to each of the prime factors.
• $X^{(1\;2)}$, the set of elements fixed by the transposition $(1\;2)$. There are $4$ possibilities for the powers of $2$: $(2^3,2^3,1)$, $(2^2,2^2,2^2)$, $(2^1,2^1,2^4)$, and $(1,1,2^6)$. There are $2$ possibilities for the powers of $3$ and $5$. So $|X^{(1\;2)}| = 16$.
• Similarly, $|X^{(1\;3)}|=|X^{(2\;3)}|=16$.
• $X^{(1\;2\;3)}$, the set of elements fixed by the $3$-cycle $(1\;2\;3)$. This means $x=y=z$ in the triple $(x,y,z)$, so there is only one such element: $(60,60,60)$. So $|X^{(1\;2\;3)}| = 1$.
• Similarly, $|X^{(1\;3\;2)}| = 1$.

So by Burnside's lemma, $$ |X/G| = \frac{2800 + 16 + 16 + 16 + 1 + 1}{6} = 475. $$

This approach is easy to generalize to counting factorizations of $xyz=n$. It generalizes to factorization with more factors as well, but then the group action is more complicated, so there are more cases to deal with, and they're individually harder to count.