Number of necklaces of 16 beads with 8 red beads, 4 green beads and 4 yellow beads The question is as stated in the title up to symmetries of $D_{16}$. I know this has to do with the following two formulas:
If $G=D_{16}$ is the group acting on the set $S$ of different necklaces, then 
1) $s\in S$, $|Orb(s)|=\frac{|G|}{|Stab(s)|}$
2) Cauchy-Frobenius Formula, so {The number of orbits of G in S}$=\frac{1}{|G|}\sum_{g\in G}f(g)$, where $f(g)=|${$s\in S : g*s=s$}$|$
I find this question is a difficult one to wrap my head around.  I'm not even quite certain what I am necessarily trying to find.  Am I trying to find the number of orbits? (I would guess this).  How do I find all the $f(g)$'s then?  I figure for example $f(1)=16$, but then if I have $f(x)$ where x is a rotation I could think of many scenarios where this differs from how you ordered it.
Any help or hints is appreciated, thanks
 A: In answering this  problem I will assume that we are  allowed to use a
computer algebra system since it is somewhat computation-intensive but
not much is gained by doing these computations with pen and paper.

Observe  that the convention  at the  OEIS is  that the  term necklace
refers to  the slots  being arranged around  a circle  with rotational
symmetry  acting  on  them.  Similarly  the  term  bracelet  indicates
reflections acting on the slots  in addition to rotations. This is the
difference  between  the cyclic  group  $C_n$  acting  on $n$  elements
(necklace) and the dihedral group $D_n$ on $n$ elements (bracelet).
In your case we have a bracelet of sixteen beads with cycle index
$$Z(D_{16}) =
1/32\,{a_{{1}}}^{16}+{\frac {9\,{a_{{2}}}^{8}}{32}}+1/16\,{a_{
{4}}}^{4}+1/8\,{a_{{8}}}^{2}+1/4\,a_{{16}}+1/4\,{a_{{1}}}^{2}{
a_{{2}}}^{7}.$$
We are interested in $$[R^8 G^4 Y^4] Z(D_{16})(R+G+Y).$$
The substituted cycle index works out to
$$1/32\, \left( R+G+Y \right) ^{16}+{\frac {9\, \left( {G}^{2}+{
R}^{2}+{Y}^{2} \right) ^{8}}{32}}+1/16\, \left( {G}^{4}+{R}^{4
}+{Y}^{4} \right) ^{4}\\+1/8\, \left( {G}^{8}+{R}^{8}+{Y}^{8}
 \right) ^{2}+1/4\,{G}^{16}+1/4\,{R}^{16}+1/4\,{Y}^{16}\\+1/4\,
 \left( R+G+Y \right) ^{2} \left( {G}^{2}+{R}^{2}+{Y}^{2}
 \right) ^{7}.$$
Here is an excerpt from the expansion of the above:
$$\ldots, 147\,{G}^{5}{Y}^{11}+72\,{G}^{4}{R}^{12}+
693\,{G}^{4}{R}^{11}Y+3843\,{G}^{4}{R}^{10}{Y}^{2}\\+12565\,{G}^
{4}{R}^{9}{Y}^{3}+28377\,{G}^{4}{R}^{8}{Y}^{4}+45150\,{G}^{4}{
R}^{7}{Y}^{5}+52850\,{G}^{4}{R}^{6}{Y}^{6},\ldots$$
The coefficient on $R^8 G^4 Y^4$ is
$$28377.$$
The Maple code for this was as follows:

with(numtheory);

pet_cycleind_cyclic :=
proc(n)
local d, s;

    s := 0;
    for d in divisors(n) do
        s := s + phi(d)*a[d]^(n/d);
    od;

    s/n;
end;

pet_cycleind_dihedral :=
proc(n)
local s;

    s := 1/2*pet_cycleind_cyclic(n);

    if(type(n, odd)) then
        s := s + 1/2*a[1]*a[2]^((n-1)/2);
    else
        s := s + 1/4*(a[1]^2*a[2]^((n-2)/2) + a[2]^(n/2));
    fi;

    s;
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;

