First we construct the generating function of the edge colorings before the symmetric group acts on the colors. This requires the cycle index of the edge permutation group $E$ of the cube. We enumerate the permutations in this group by their cycle structure.
First, there is the identity, which contributes $$a_1^{12}.$$
There are rotations about an axis (three of them) passing through the center of two opposite faces, which gives $$3\times (a_2^6 + 2 a_4^3).$$
Rotations about a diagonal passing through opposite vertices give
$$4\times 2 a_3^4.$$
Finally there are 180 degree rotations about an axis passing through the centers of opposite edges, which gives $$6\times a_1^2 a_2^5.$$
This gives the cycle index
$$Z(E) = \frac{1}{24}
\left(a_1^{12} + 3 a_2^6 + 6 a_4^3 + 8 a_3^4 + 6 a_1^2 a_2^5\right).$$
Substituting into $Z(E)$ we obtain the following formula for the number of $N$-colorings of the edges of a cube the value
$$\frac{1}{24}(N^{12} + 3 N^6 + 6 N^3 + 8 N^4 + 6 N^7).$$
This gives the following sequence
$$1, 218, 22815, 703760, 10194250, 90775566, 576941778, 2863870080, 11769161895,\ldots$$
which is OEIS A060530.
The substituted cycle index for 3-colorings as asked for by the OP is
$$Z(E)(R+G+B)=
1/24\, \left( R+G+B \right) ^{12}+1/8\, \left( {R}^{2}+{G}^{2}+{B}^{2}
\right) ^{6}\\+1/4\, \left( {R}^{4}+{G}^{4}+{B}^{4} \right) ^{3}+1/3\, \left(
{R}^{3}+{G}^{3}+{B}^{3} \right) ^{4}+1/4\, \left( R+G+B \right) ^{2} \left( {
R}^{2}+{G}^{2}+{B}^{2} \right) ^{5}$$
which yields
$$[R^4 G^4 B^4]Z(E)(R+G+B)=1479.$$
Addendum as of Tue Dec 17 22:50:41 CET 2013. The OP also mentioned in his query that as an extra twist we may consider edge colorings of the cube where the symmetric group acts on the colors, so that two colorings are not only equivalent when there is a rotation that transforms one into the other, but also a rotation combined with a permutation of the colors from the symmetric group. The technique for this kind of problem is called Power Group Enumeration Theorem.
There are at least two possible approaches here. There is the classic approach from the book Graphical Enumeration by Frank Harary and Edgar Palmer, a landmark text that continues to be relevant today. They have a chapter on this namely chapter six, section four "Power Group Enumeration / Graphs with colored lines." I will not typeset their formulas here as I would not add anything new and the text is available at any major library. Let me just say that they treat the following scenario: enumerate under the power group the colorings of various graphs, where an edge maybe on or off and has $N$ possible colors in the on state. I implemented their formula in Maple and the result was the following sequence for edge colorings of the cube under the symmetric group acting on the colors using at most $N$ colors:
$$1, 114, 3891, 29854, 87981, 143797, 170335, 177160, 178153, 178243,\ldots$$
This is the Maple code for the Harary/Palmer computation.
with(numtheory);
with(group):
with(combinat):
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_power :=
proc(N, ZA)
local beta, res, indvars, v, t, d, m, pot, ZB;
ZB := pet_cycleind_symm(N);
if N=1 then ZB := [ZB]; fi;
indvars := indets(ZA);
res := 0;
for beta in ZB do
t := [];
for v in indvars do
pot := op(1, v);
m := 0;
for d in divisors(pot) do
m := m +
d*degree(beta, op(0,v)[d]);
od;
t := [op(t), v=1+m*z^pot];
od;
res := res+lcoeff(beta)*subs(t, ZA);
od;
res;
end;
cube_edge_ind :=
1/24*(a[1]^12+3*a[2]^6+6*a[4]^3+8*a[3]^4+6*a[1]^2*a[2]^5);
vgf :=
proc(N)
option remember;
local gf;
gf := pet_varinto_power(N, cube_edge_ind);
expand(gf);
end;
v := N -> coeff(vgf(N), z, 12);
The second approach uses a different Power Group Enumeration (PGET) formula as presented by Harald Fripertinger in his brilliant paper "Enumeration in Musical Theory" from 1993. Again I will not typeset the formulas here as I would not be adding anything new. The clinch is that this paper presents a so-called polynomial form of the PGET (in the first chapter) which makes it possible not only to calculate the total number of colorings but also to classify them according to the distribution of colors used. I have implemented this in Maple and I obtain for the sequence of colorings once more the values (that is as it ought to be):
$$1, 114, 3891, 29854, 87981, 143797, 170335,\ldots$$
Now, however, we also get the distributions of colors. For two colors we have
$$29\,{\it Q6\_Q6}+{\it Q12}+5\,{\it Q2\_Q10}+27\,{\it Q4\_Q8}+{\it Q1\_Q11}+13\,{\it
Q3\_Q9}+38\,{\it Q5\_Q7} $$
and for three colors
$$ {\it Q12}+38\,{\it Q5\_Q7}+27\,{\it Q4\_Q8}+{\it Q1\_Q11}+5\,{\it Q2\_Q10}+13\,{\it
Q3\_Q9}+370\,{\it Q2\_Q5\_Q5}\\+29\,{\it Q6\_Q6}+600\,{\it Q2\_Q4\_Q6}+236\,{\it
Q1\_Q5\_Q6}+85\,{\it Q1\_Q3\_Q8}+77\,{\it Q2\_Q2\_Q8}\\+30\,{\it Q1\_Q2\_Q9}+340\,{\it
Q2\_Q3\_Q7}+5\,{\it Q1\_Q1\_Q10}+170\,{\it Q1\_Q4\_Q7}\\+412\,{\it Q3\_Q3\_Q6}+282\,{\it
Q4\_Q4\_Q4}+1170\,{\it Q3\_Q4\_Q5}$$
In particular we are now in a position to answer the query that prompted the post which is how many colorings there are using each of the three colors four times and the answer is
$$282$$
Here are the colorings using at most four colors:
$${\it Q12}+38\,{\it Q5\_Q7}+27\,{\it Q4\_Q8}+{\it Q1\_Q11}+5\,{\it Q2\_Q10}+13\,{\it
Q3\_Q9}+370\,{\it Q2\_Q5\_Q5}\\+698\,{\it Q3\_Q3\_Q3\_Q3}+29\,{\it Q6\_Q6}+3480\,{\it
Q1\_Q2\_Q4\_Q5}+3510\,{\it Q2\_Q2\_Q3\_Q5}\\+600\,{\it Q2\_Q4\_Q6}+2330\,{\it
Q1\_Q3\_Q3\_Q5}+2915\,{\it Q1\_Q3\_Q4\_Q4}+600\,{\it Q1\_Q1\_Q4\_Q6}\\+510\,{\it
Q1\_Q2\_Q2\_Q7}+2266\,{\it Q2\_Q2\_Q4\_Q4}+236\,{\it Q1\_Q5\_Q6}+85\,{\it Q1\_Q3\_Q8}\\+
77\,{\it Q2\_Q2\_Q8}+30\,{\it Q1\_Q2\_Q9}+2320\,{\it Q1\_Q2\_Q3\_Q6}+340\,{\it
Q1\_Q1\_Q3\_Q7}\\+13\,{\it Q1\_Q1\_Q1\_Q9}+626\,{\it Q2\_Q2\_Q2\_Q6}+340\,{\it Q2\_Q3\_Q7
}+5\,{\it Q1\_Q1\_Q10}\\+5850\,{\it Q2\_Q3\_Q3\_Q4}+135\,{\it Q1\_Q1\_Q2\_Q8}+170\,{\it
Q1\_Q4\_Q7}+412\,{\it Q3\_Q3\_Q6}\\+282\,{\it Q4\_Q4\_Q4}+370\,{\it Q1\_Q1\_Q5\_Q5}+1170
\,{\it Q3\_Q4\_Q5} $$
The code for the Fripertinger computation follows.
with(numtheory);
with(group):
with(combinat):
pet_disjcyc :=
proc(p)
local dc, pos;
dc := convert(p, 'disjcyc');
for pos to nops(p) do
if p[pos] = pos then
dc := [op(dc), [pos]];
fi;
od;
dc;
end;
pet_varinto_power :=
proc(N, ZA)
local indvars, delta, v, pot, s, t, poly, k, degs,
res, cyc;
indvars := indets(ZA);
res := 0;
for delta in permute(N) do
t := [];
for v in indvars do
pot := op(1, v);
s := 0;
for cyc in pet_disjcyc(delta) do
if pot mod nops(cyc) = 0 then
poly :=
mul(cat(`X`, cyc[k])^(pot/nops(cyc)),
k=1..nops(cyc));
s := s + nops(cyc)*poly;
fi;
od;
t := [op(t), v=s];
od;
res := res+subs(t, ZA);
od;
degs :=
proc(t) local q;
local vs, v, dl, d, sym;
vs := indets(t);
dl := [];
for v in vs do
dl := [op(dl), degree(t, v)];
od;
sym := [];
for d in sort(dl) do
if nops(sym) = 0 then
sym := cat(`Q`, d);
else
sym := cat(sym, `_`, `Q`, d);
fi;
od;
return lcoeff(t)*sym;
end;
if N=1 then
return cat(`Q`, degree(res, `X1`));
fi;
map(degs, expand(res))/N!;
end;
cube_edge_ind :=
1/24*(a[1]^12+3*a[2]^6+6*a[4]^3+8*a[3]^4+6*a[1]^2*a[2]^5);
v :=
proc(N)
option remember;
local dist;
dist := pet_varinto_power(N, cube_edge_ind);
map(lcoeff, dist);
end;
Thanks go to the OP for asking such a fascinating question. Here is a link to a chain of Polya Counting Computations.
Addendum as of Sat Dec 21 03:20:27 CET 2013. The above algorithm can also be used on generic graphs rather than just cubes. With two colors, we get the sequence
$$1, 2, 6, 18, 78, 522, 6178, 137352, 6002584, 509498932, 82545586656, 25251015686776,\ldots$$
which is OEIS A007869.
Three interchangeable colors produces the sequence
$$1, 3, 15, 142, 4300, 384199, 98654374, 70130880569, \\136638863494089, 730439999032117301,\ldots$$
which does not yet have an OEIS entry.
With four interchangeable colors we get
$$1, 3, 22, 513, 67685, 37205801, 74992370359, 543437207831908, 14224090440652751128,\ldots $$
The computation of the relevant cycle indices is show at this MSE link. This MSE link II shows Power Group Enumeration with the flip group acting on the slots.
Addendum Dec 13 2016. The reader may be interested to know that
these generating functions can also be computed from first principles
namely the Burnside lemma as was done at the following MSE
link. The algorithm
is documented there and in the linked-to pages. I present this
algorithm here for reference. The output up to sort order is the same
as from the Fripertinger computation.
with(combinat);
pet_flatten_term :=
proc(varp)
local terml, d, cf, v;
terml := [];
cf := varp;
for v in indets(varp) do
d := degree(varp, v);
terml := [op(terml), seq(v, k=1..d)];
cf := cf/v^d;
od;
[cf, terml];
end;
pet_autom2cyclesA :=
proc(src, aut)
local numa, numsubs;
local marks, pos, data, item, cpos, clen;
numsubs := [seq(src[k]=k, k=1..nops(src))];
numa := subs(numsubs, aut);
marks := Array([seq(true, pos=1..nops(aut))]);
pos := 1; data := table();
while pos <= nops(aut) do
if marks[pos] then
clen := 0; item := []; cpos := pos;
while marks[cpos] do
marks[cpos] := false;
item := [op(item), aut[cpos]];
clen := clen+1;
cpos := numa[cpos];
od;
if type(data[clen], 'list') then
data[clen] :=
[op(data[clen]), item];
else
data[clen] := [item];
fi;
fi;
pos := pos+1;
od;
return data;
end;
pet_flatten_automcycs :=
proc(tbl)
local res, sz, cyc;
res := [];
for sz in [indices(tbl, 'nolist')] do
for cyc in tbl[sz] do
res := [op(res), cyc];
od;
od;
res;
end;
cube_edge_cind :=
1/24*(a[1]^12 + 8*a[3]^4 + 6*a[4]^3 + 3*a[2]^6 + 6*a[1]^2*a[2]^5);
cube_edge_colorings :=
proc(n)
option remember;
local idx_colors, res, a, b,
flat_a, flat_b, cyc_a, cyc_b, len_a, len_b, p, q,
src, perm, acycs, gf, msetgf, term, deg, rep;
res := 0; src := [seq(el, el=1..n)];
for a in cube_edge_cind do
flat_a := pet_flatten_term(a);
perm := firstperm(n);
while type(perm, `list`) do
acycs := pet_autom2cyclesA(src, perm);
flat_b := pet_flatten_automcycs(acycs);
p := 1;
for cyc_a in flat_a[2] do
len_a := op(1, cyc_a);
q := 0;
for cyc_b in flat_b do
len_b := nops(cyc_b);
if len_a mod len_b = 0 then
q := q +
len_b*mul(Q[el], el in cyc_b)
^(len_a/len_b);
fi;
od;
p := p*q;
od;
res := res + p*flat_a[1];
perm := nextperm(perm);
od;
od;
gf := expand(res/n!);
msetgf := 0;
for term in gf do
rep := 0;
for deg in
sort(map(Q-> degree(term, Q),
convert(indets(term), `list`)))
do
if type(rep, `integer`) then
rep := cat(`Q`, deg);
else
rep := cat(rep, `_Q`, deg);
fi;
od;
msetgf := msetgf + lcoeff(term)*rep;
od;
msetgf;
end;
