Consider the problem of $4n$ balls with $n$ balls of each of the four
colors being distributed into four indistinguishable baskets where
each basket holds exactly $n$ balls. The naive approach here would be
to use the Polya Enumeration Theorem (twice). Surprisingly enough this
is sufficient to compute the initial segment of the sequence using the
recurrence by Lovasz for the cycle index $Z(S_n)$ of the multiset
operator $\def\textsc#1{\dosc#1\csod} \def\dosc#1#2\csod{{\rm
#1{\small #2}}} \textsc{MSET}_{=n}$ on $n$ slots, 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.$$
This recurrence lets us calculate the cycle index $Z(S_n)$ very
easily. The answer is then given by
$$[A^n B^n C^n D^n] Z(S_4)(Z(S_n)(A+B+C+D)).$$
Using Maple and a reasonable amount of computational resources this
yields the sequence
$$1, 17, 93, 465, 1746, 5741, 16238, 41650, 97407, 212412, 434767,
\\ 845366, 1569344, 2801696, 4828140, 8069053,
\\ 13114785, 20796651, 32242621, 48986553, 73052382, 107114645,
\\ 154621230, 220021932, 308940815,\ldots$$
In particular the value for $n=10$ is given by
$$212412.$$
This is OEIS A253259 where we discover a
variation of the problem definition that confirms the validity of
these results. Oddly enough no recurrence relation or other indication
of how these numbers were computed is given in the OEIS entry. Perhaps
we will see a recurrence now that there are enough test data to verify
its correctness, if indeed it exists.
The Maple code for the above is quite straightforward.
with(combinat);
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;
pet_cycleind_symm :=
proc(n)
local l;
option remember;
if n=0 then return 1; fi;
expand(1/n*add(a[l]*pet_cycleind_symm(n-l), l=1..n));
end;
V :=
proc(n)
option remember;
local comb, gf, var;
comb := pet_varinto_cind(A+B+C+D, pet_cycleind_symm(n));
gf :=
expand(pet_varinto_cind(comb, pet_cycleind_symm(4)));
for var in [A,B,C,D] do
gf := coeff(gf, var, n);
od;
gf;
end;
Addendum. I realized I had overlooked an additional OEIS link to
A257463 when I wrote the above several
hours ago. It proposes a simple algorithm which uses the isomorphism
between the problem and factorizations of $p_1^n p_2^n p_3^n p_4^n$
into four factors $q$ all of which have $\Omega(q) = n.$ The algorithm
generates all of these using the observation that uniqueness of
factorizations can be guaranteed by generating the factors in
non-increasing order. It uses memoization to speed this up.
Nonetheless when I pasted it verbatim into Maple and tried it on
several test cases it performed very poorly compared to what I have
above. I will therefore keep the post for now. The quest continues.
Addendum II. As per request by @WillOrrick I am posting the code
for the general problem of $k$ colors.
V :=
proc(n, k)
option remember;
local base, comb, gf, var;
base := add(Q[p], p=1..k);
comb := pet_varinto_cind(base, pet_cycleind_symm(n));
gf :=
expand(pet_varinto_cind(comb, pet_cycleind_symm(k)));
for var in [seq(Q[p], p=1..k)] do
gf := coeff(gf, var, n);
od;
gf;
end;
We thus obtain for five colors the sequence
$$1, 73, 1417, 19834, 190131, 1398547, 8246011, 40837569, 174901563,
\\ 664006236, 2274999093, 7139338769, 20758868781, 56466073587,
\\ 144806582536, 352420554194, 818441723112, 1822255658908,\ldots$$
Addendum III. Granting Maple several hours of computation time and
5GB of memory we get for four colors:
$$1, 17, 93, 465, 1746, 5741, 16238, 41650, 97407, 212412, 434767,
\\ 845366, 1569344, n2801696, 4828140, 8069053, 13114785, 20796651,
\\ 32242621, 48986553, 73052382, 107114645, 154621230, 220021932,
\\ 308940815, 428492880, 587520315, 797019526, 1070458096, 1424339518,
\\ 1878618620, 2457435561, 3189651885, 4109787687, 5258703597,\ldots$$
This confirms the generating function by @WillOrrick.
