Cube stack problem

We will solve this problem with the Polya Enumeration Theorem. Evidently this is equivalent to counting two-colorings of the small cubes (leaving out a cube is just like coloring it a second, different color). This requires the cycle index of the permutation group $G$ of the small cubes inside the large one. It is important to note that we restrict ourselves to rotations. If we had a graph instead of a real-world object there would be mirror reflections to consider as well (reflections about a plane passing through the center of the cube and parallel to two opposite sides). This cannot be done with a physical concrete object. We now enumerate the permutations of $Z(G)$ by their cycle structure, considering the types of rotations in turn.

Start with the identity, which contributes $$a_1^{27}.$$

Now there are eight rotations about axes passing through opposite corners of the cube (two rotations per axis, of which there are four), which fix the three cubes on the diagonal and which contribute $$8\times a_1^3 a_3^8.$$ There are three rotations about axes passing through the centers of opposite faces, which again fix the cubes on said axis and which give $$3\times (2 a_1^3 a_4^6 + a_1^3 a_2^{12}).$$ Finally we have six rotations about axes passing through the centers of opposite edges, which fix the cubes on that axis and give $$6\times a_1^3 a_2^{12}.$$ This gives the cycle index $$Z(G) = \frac{1}{24} \left(a_1^{27} + 8\times a_1^3 a_3^8 + 3\times (2 a_1^3 a_4^6 + a_1^3 a_2^{12}) + 6\times a_1^3 a_2^{12}\right)$$ which is $$Z(G) = \frac{1}{24} \left(a_1^{27} + 8\times a_1^3 a_3^8 + 6\times a_1^3 a_4^6 + 9\times a_1^3 a_2^{12}\right).$$ Now this gives $$Z(G)(1+z) = 1/24\, \left( 1+z \right) ^{27}+1/3\, \left( 1+z \right) ^{3} \left( 1+{z}^{3} \right) ^{8}\\+1/4\, \left( 1+z \right) ^{3} \left( 1+{z}^{4} \right) ^{6}+3/8\, \left( 1+z \right) ^{3} \left( 1+{z}^{2} \right) ^{12}$$ which is $$Z(G)(1+z) = {z}^{27}+4\,{z}^{26}+22\,{z}^{25}+139\,{z}^{24}+779\,{z}^{23}+3455\,{z}^{22}\\+12507\, {z}^{21}+37303\,{z}^{20}+92968\,{z}^{19}+195963\,{z}^{18}+352433\,{z}^{17}+544382\,{ z}^{16}\\+725612\,{z}^{15}+837184\,{z}^{14}+837184\,{z}^{13}+725612\,{z}^{12}+544382\, {z}^{11}+352433\,{z}^{10}\\+195963\,{z}^{9}+92968\,{z}^{8}+37303\,{z}^{7}+12507\,{z}^{ 6}+3455\,{z}^{5}\\+779\,{z}^{4}+139\,{z}^{3}+22\,{z}^{2}+4\,z+1,$$ so the answers are four, twenty-two, and one hundred thirty nine respectively.

This MSE link computes the full face permutation group of the n-dimensional hypercube.

For the sake of completeness here is a companion answer that includes reflections in addition to rotations. This was result was not computed manually and the reader is invited to give a combinatorial proof for the terms of the cycle index of the full symmetry of the cube. (This should be doable using the fact that the additional permutations compared to rotations only are obtainable from the former by performing an inversion, but it may need building a paper model.)

Here we used the known fact that the symmetries of the hypercube on $2^n$ vertices are combinations of bit flips and bit permutations, where the vertices are labeled with $n$ bit strings and two vertices are adjacent if they differ in one bit. A bit flip flips some fixed subset of the $n$ bits and a permutation permutes them. Each combination of these two yields a vertex permutation that is guaranteed to be an automorphism. This makes it easy to compute the cycle index because we can enumerate the $2^n \times n!$ permutations relatively quickly for small $n$ (as opposed to searching the space of all $(2^n)!$ permutations). Note that hypercubes of higher dimensions are definitely tractable with this method and the reader is urged to give this a try.

The cycle index $Z(H)$ for the full symmetries of the $3\times 3\times 3$ cube was found to be $$Z(H) = \frac{1}{48} \left( {a_{{1}}}^{27}+9\,{a_{{1}}}^{9}{a_{{2}}}^{9}+9\,{a_{{1}}}^{3}{a_{{2}}}^{12}+a_{{1}}{a_{{2}}}^{13} +8\,{a_{{1}}}^{3}{a_{{3}}}^{8}+6\,{a_{{1}}}^{3}{a_{{4}}}^{6}\\+6\,a_{{1}}a_{{2}}{a_{{4}}}^{6}+8\,a_ {{1}}a_{{2}}{a_{{6}}}^{4}\right).$$

The generating functions for two-color patterns or alternatively, removed cubes, turns out to be the following $$Z(H)(1+z) = 1/48\, \left( 1+z \right) ^{27}+3/16\, \left( 1+z \right) ^{9} \left( 1+{z}^{2} \right) ^{9}+1/6 \, \left( 1+z \right) ^{3} \left( 1+{z}^{3} \right) ^{8}\\+3/16\, \left( 1+z \right) ^{3} \left( 1+ {z}^{2} \right) ^{12}+1/8\, \left( 1+z \right) ^{3} \left( 1+{z}^{4} \right) ^{6}+1/6\, \left( 1+ {z}^{6} \right) ^{4} \left( 1+{z}^{2} \right) \left( 1+z \right)\\ +1/8\, \left( 1+{z}^{4} \right) ^{6} \left( 1+{z}^{2} \right) \left( 1+z \right) +1/48\, \left( 1+z \right) \left( 1+{ z}^{2} \right) ^{13}.$$

Expansion gives $$Z(H)(1+z) = {z}^{27}+4\,{z}^{26}+20\,{z}^{25}+101\,{z}^{24} +483\,{z}^{23}+1956\,{z}^{22}+6748\,{z}^{21}\\+19587 \,{z}^{20}+48086\,{z}^{19}+100446\,{z}^{18}+179686\,{z}^{17}+276646\,{z}^{16}+368072\,{z}^{15}\\+ 424308\,{z}^{14}+424308\,{z}^{13}+368072\,{z}^{12}+276646\,{z}^{11} +179686\,{z}^{10}+100446\,{z}^ {9}\\+48086\,{z}^{8}+19587\,{z}^{7}+6748\,{z}^{6}+1956\,{z}^{5} +483\,{z}^{4}+101\,{z}^{3}+20\,{z}^{ 2}+4\,z+1.$$ The sequence of unique colorings of the $3\times 3\times 3$ with $n$ colors starts out like this: $$1, 2852288, 158942078604, 375313061509120, 155221150215953125, 21322735154793366336,$$ which does not yet have an OEIS entry.

This is the Maple code that was used in the computation. Questions and comments are welcome as there usually is room for imrovement in the first version of a computer program. I could have hard-coded the representations of the small cubes as opposed to calculating them but that would not generalize to hypercubes of higher dimensions.

with(numtheory);
with(group):
with(combinat):

pet_autom2cycles :=
proc(src, aut)
        local numa, numsubs;
        local marks, pos, cycs, cpos, clen;

        numsubs := [seq(src[k]=k, k=1..nops(src))];
        numa := subs(numsubs, aut);

        marks := [seq(true, pos=1..nops(aut))];

        cycs := []; pos := 1;

        while pos <= nops(aut) do
              if marks[pos] then
                 clen := 0; cpos := pos;

                 while marks[cpos] do
                       marks[cpos] := false;
                       cpos := numa[cpos];
                       clen := clen+1;
                 od;

                 cycs := [op(cycs), clen];
              fi;

              pos := pos+1;
        od;

        return mul(a[cycs[k]], k=1..nops(cycs));
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;

cube_vsyms :=
proc()
         option remember;
         local aut, adj, p, bitstr, bits, b, vsyms, vsym, flip,
         doflip, doperm, bits2ind, v1, v2, idx, idxterm,
         corners, edgemids, centers, cen, censet, c;

         bitstr := [];
         for b from 0 to 7 do
             bits := convert(8+b, base, 2);
             bitstr := [op(bitstr), [seq(bits[k], k=1..3)]];
         od;

         bits2ind := b -> 1+b[1]+2*b[2]+4*b[3];

         vsyms := [];

         for flip from 0 to 7 do
             bits := convert(8+flip, base, 2);

             doflip :=
             proc(v)
                 local res, bpos;
                 res := [];

                 for bpos to 3 do
                     if bits[bpos] = 1 then
                         res := [op(res), 1-v[bpos]];
                     else
                         res := [op(res), v[bpos]];
                     fi;
                 od;
             end;

             for p in permute(3) do
                 doperm :=
                 proc(v)
                     return [v[p[1]], v[p[2]], v[p[3]]];
                 end;

                 vsym := map(doflip, bitstr);
                 vsym := map(doperm, vsym);

                 vsyms := [op(vsyms), map(bits2ind, vsym)];
             od;
         od;

         corners := [seq(k,k=1..8)];

         edgemids := [];
         for v1 to 8 do
             for v2 from v1+1 to 8 do
                 adj :=
                 proc(b1, b2)
                     local df;

                     df := 0;
                     if b1[1] <> b2[1] then df := df+1 fi;
                     if b1[2] <> b2[2] then df := df+1 fi;
                     if b1[3] <> b2[3] then df := df+1 fi;

                     return df
                 end;

                 if adj(bitstr[v1], bitstr[v2]) = 1 then
                     edgemids := [op(edgemids), {v1, v2}];
                 fi;
             od;
         od;

         centers := [];
         cen := [[0,0], [0,1], [1, 0], [1, 1]];
         for b to 3 do
             censet := [];
             for c to 4 do
                 censet :=
                 [op(censet),
                  [seq(cen[c][k], k=1..b-1),
                   0,
                   seq(cen[c][k], k=b..2)]];
             od;
             centers :=
             [op(centers),
              convert(map(bits2ind, censet), set)];

             censet := [];
             for c to 4 do
                 censet :=
                 [op(censet),
                  [seq(cen[c][k], k=1..b-1),
                   1,
                   seq(cen[c][k], k=b..2)]];
             od;
             centers :=
             [op(centers),
              convert(map(bits2ind, censet), set)];
         od;

         idx := 0;

         for aut in vsyms do
             idxterm := 1;

             p := [seq(k=aut[k], k=1..nops(aut))];

             idxterm :=
             idxterm *
             pet_autom2cycles(corners, subs(p, corners));

             idxterm :=
             idxterm *
             pet_autom2cycles(edgemids, subs(p, edgemids));

             idxterm :=
             idxterm *
             pet_autom2cycles(centers, subs(p, centers));

             idx := idx+idxterm*a[1];
         od;

         idx/nops(vsyms);
end;


v :=
proc(n)
         option remember;
         local gf, rep;

         rep := add(cat('X', k), k=1..n);
         gf := expand(pet_varinto_cind(rep, cube_vsyms()));

         subs([seq(cat('X',k)=1, k=1..n)], gf);
end;

For those who are interested here is a version of the Maple program that has hard-coded vertices, edges and faces. It is not portable to higher dimension hypercubes but as the code is reduced considerably and readability profits, with the focus being on the essential.

with(numtheory);
with(group):
with(combinat):

pet_autom2cycles :=
proc(src, aut)
        local numa, numsubs;
        local marks, pos, cycs, cpos, clen;

        numsubs := [seq(src[k]=k, k=1..nops(src))];
        numa := subs(numsubs, aut);

        marks := [seq(true, pos=1..nops(aut))];

        cycs := []; pos := 1;

        while pos <= nops(aut) do
              if marks[pos] then
                 clen := 0; cpos := pos;

                 while marks[cpos] do
                       marks[cpos] := false;
                       cpos := numa[cpos];
                       clen := clen+1;
                 od;

                 cycs := [op(cycs), clen];
              fi;

              pos := pos+1;
        od;

        return mul(a[cycs[k]], k=1..nops(cycs));
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;

cube_vsyms :=
proc()
         option remember;
         local aut, adj, p, bitstr, bits, b, vsyms, vsym, flip,
         doflip, doperm, bits2ind, v1, v2, idx, idxterm,
         corners, edgemids, centers;

         bitstr := [];
         for b from 0 to 7 do
             bits := convert(8+b, base, 2);
             bitstr := [op(bitstr), [seq(bits[k], k=1..3)]];
         od;

         bits2ind := b -> 1+b[1]+2*b[2]+4*b[3];

         vsyms := [];

         for flip from 0 to 7 do
             bits := convert(8+flip, base, 2);

             doflip :=
             proc(v)
                 local res, bpos;
                 res := [];

                 for bpos to 3 do
                     if bits[bpos] = 1 then
                         res := [op(res), 1-v[bpos]];
                     else
                         res := [op(res), v[bpos]];
                     fi;
                 od;
             end;

             for p in permute(3) do
                 doperm :=
                 proc(v)
                     return [v[p[1]], v[p[2]], v[p[3]]];
                 end;

                 vsym := map(doflip, bitstr);
                 vsym := map(doperm, vsym);

                 vsyms := [op(vsyms), map(bits2ind, vsym)];
             od;
         od;

         corners := [1, 2, 3, 4, 5, 6, 7, 8];
         edgemids := [{1, 2}, {1, 3}, {1, 5}, {2, 4}, {2, 6},
                      {3, 4}, {3, 7}, {4, 8}, {5, 6}, {5, 7},
                      {6, 8}, {7, 8}];
         centers := [{1, 3, 5, 7}, {2, 4, 6, 8}, {1, 2, 5, 6},
                     {3, 4, 7, 8}, {1, 2, 3, 4}, {5, 6, 7, 8}];

         idx := 0;

         for aut in vsyms do
             idxterm := 1;

             p := [seq(k=aut[k], k=1..nops(aut))];

             idxterm :=
             idxterm *
             pet_autom2cycles(corners, subs(p, corners));

             idxterm :=
             idxterm *
             pet_autom2cycles(edgemids, subs(p, edgemids));

             idxterm :=
             idxterm *
             pet_autom2cycles(centers, subs(p, centers));

             idx := idx+idxterm*a[1];
         od;

         idx/nops(vsyms);
end;


v :=
proc(n)
         option remember;
         local gf, rep;

         rep := add(cat('X', k), k=1..n);
         gf := expand(pet_varinto_cind(rep, cube_vsyms()));

         subs([seq(cat('X',k)=1, k=1..n)], gf);
end;

Here is another interesting MSE cycle index computation I. This MSE cycle index computation II could be qualified as rather exotic.