This problem is a straightforward application of the Polya Enumeration
Theorem. Suppose we treat the case of $r$ red balls, $g$ green balls
and $b$ blue balls and $n$ indistinguishable bins where no bins are
left empty.
Recall the recurrence by Lovasz for the cycle index $Z(P_n)$ of the
set operator $\def\textsc#1{\dosc#1\csod} \def\dosc#1#2\csod{{\rm
#1{\small #2}}}\textsc{SET}_{=n}$ on $n$ slots, which is $$Z(P_n) =
\frac{1}{n} \sum_{l=1}^n (-1)^{l-1} a_l Z(P_{n-l})
\quad\text{where}\quad Z(P_0) = 1.$$
We need to employ the set operator because the elements of a distribution
are supposed to be unique.
We have for example,
$$Z(P_3) = 1/6\,{a_{{1}}}^{3}-1/2\,a_{{2}}a_{{1}}+1/3\,a_{{3}}$$
and
$$Z(P_4) = 1/24\,{a_{{1}}}^{4}-1/4\,a_{{2}}{a_{{1}}}^{2}
+1/3\,a_{{3}}a_{{1}}+1/8\,{a_{{2}}}^{2}-1/4\,a_{{4}}.$$
Applying PET it now follows almost by inspection that the
desired count is given by using the repertoire
$$-1 + \sum_{q=0}^r R^q \sum_{q=0}^g G^q
\sum_{q=0}^b B^q$$
with $Z(P_n)$ to get
$$[R^r G^g B^b] Z\left(P_n; -1 +
\sum_{q=0}^r R^q \sum_{q=0}^g G^q \sum_{q=0}^b B^q \right).$$
The minus one term cancels empty bins.
The answer for eight red, six green and seven blue balls in four
indistinguishable bins turns out to be
$$60040.$$
The following Maple code implements two routines, q1 and
q2. The first of these computes the value for $(r,g,b)$ and $n$ by
brute force (enumerate all configurations) and can be used to verify
the correctness of the PET formula for small values of the parameters.
The second one uses the PET for instant computation of the desired
coefficient.
Observe that we can compute values that are utterly out of reach of
brute force enumeration, e.g. with ten balls of each color and five
bins we obtain
$$7098688.$$
With twenty balls of each color and six bins we get
$$194589338219.$$
with(combinat);
pet_cycleind_set :=
proc(n)
option remember;
if n=0 then return 1; fi;
expand(1/n*
add((-1)^(l-1)*a[l]*pet_cycleind_set(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;
allparts :=
proc(val, size)
option remember;
local res, els, p, pp, q;
res := [];
for p in partition(val) do
if nops(p) <= size then
pp := [op(p), 0$(size-nops(p))];
res := [op(res), pp];
fi;
od;
res;
end;
q1 :=
proc(RC, GC, BC, n)
option remember;
local p, q, pr, pg, pb, res,
sr, sg, sb, dist;
pr := [];
for p in allparts(RC, n) do
pr := [op(pr), op(permute(p))];
od;
pg := [];
for p in allparts(GC, n) do
pg := [op(pg), op(permute(p))];
od;
pb := [];
for p in allparts(BC, n) do
pb := [op(pb), op(permute(p))];
od;
res := {};
for sr in pr do
for sg in pg do
for sb in pb do
dist :=
{seq(R^sr[pos]*G^sg[pos]*B^sb[pos],
pos=1..n)};
if nops(dist) = n and
not member(1, dist) then
res := res union {dist};
fi;
od;
od;
od;
nops(res);
end;
q2 :=
proc(RC, GC, BC, n)
option remember;
local q, rep, sind;
rep := add(R^q, q=0..RC)
* add(G^q, q=0..GC)
* add(B^q, q=0..BC);
sind := pet_varinto_cind(rep-1,
pet_cycleind_set(n));
coeff(coeff(coeff(expand(sind), R, RC),
G, GC), B, BC);
end;
Applying $S_2$ here will be quite convoluted.
As explained in another question of yours,
The number of ways of getting exactly k different values when you throw n fair r-sided dice is
$\binom{r}{k}\cdot S_2(n,k)\cdot k! $
so to get exactly $6$ faces in $12$ throws of a normal die, $\binom66\cdot S_2(12,6)\cdot6!$ ways,
but this would include $11$ patterns, viz.
$\boxed 7\boxed 1\boxed 1\boxed 1\boxed 1\boxed 1\;\;\boxed 6\boxed 2\boxed 1\boxed 1\boxed 1\boxed 1\;\;\boxed 5\boxed 3\boxed 1\boxed 1\boxed 1\boxed 1\;\; \boxed 5\boxed 2\boxed 2\boxed 1\boxed 1\boxed 1\;\;\boxed 4\boxed 4\boxed 1\boxed 1\boxed 1\boxed 1$
$\boxed 4\boxed 3\boxed 2\boxed 1\boxed 1\boxed 1\;\;\boxed 4\boxed 2\boxed 2\boxed 2\boxed 1\boxed 1\;\;\boxed 3\boxed 3\boxed 3\boxed 1\boxed 1\boxed 1 \;\;\boxed 3\boxed 3\boxed 2\boxed 2\boxed 1\boxed 1\;\;\boxed 3\boxed 2\boxed 2\boxed 2\boxed 2\boxed 1$
$\boxed 2\boxed 2 \boxed 2 \boxed 2 \boxed 2 \boxed 2$
and we are interested only in the last pattern.
Now the number of ways for placing $12$ distinguishable balls in $6$ bins, the bins being indistinguishable except by occupancy, the multinomial coefficient has to be divided by ${p! q!..}$ where $p$ bins each have $p_n$ balls, $q$ bins each have $q_n$ balls, etc, e.g. for $\boxed 4\boxed 2\boxed 2\boxed 2\boxed 1\boxed 1$, there are $\binom{12}{4,2,2,2,1,1}/(3!2!)$ ways.
$S_2(12,6)= 1,323,652$ From this, if we subtract number of ways for $10$ unwanted patterns, we get
$1323652 - \left[\binom{12}{7,1,1,1,1,1}/ 5!+\binom{12}{6,2,1,1,1,1,}/4! ...+\binom{12}{5,2,2,1,1,1}/(2!3!) + ... \binom{12}{3,2,2,2,2,1}/4!\right] = 10,395$
$10,395$ is the number of ways to place$2$ each of $12$ distinguishable balls in $6$ indistinguishable bins,
So $10,395\times 6! = 7,484,400$, the desired answer.
But this looks like trying to catch your nose by bringing your hand from behind your neck !
Best Answer
As far as i have understood the boxes are distinguishable and we want to disperse balls , lets call them $(a,a), (b,b),(c,c),(d,d),(e,e),(f,f),(g,g),(h,h)$ without restriction.
In this question ,dispersing different ball pairs are independent from each other , so disperse them separately and multiply them such that
Dispersing $(a,a)$ to three distingusihable bins : $$C(3+2-1,2)=6$$ ways
Dispersing $(b,b)$ to three distingusihable bins: $$C(2+3-1,2)=6$$ ways
This process goes to up to $(h,h)$
Then , $$6 \times 6\times 6 \times 6\times 6 \times 6\times 6\times6 = 6^8=1679616$$