NB. Very inefficent wrt memory and somewhat inefficent wrt speed,
NB. but no need to optimize on a modern machine.

NB. a cube is a list of 6 sides (each atoms): 
NB. left front right back top bottom.
NB. the first four of these sides are visible on the sides in a tower.
NB. a towers is a list of 4 cubes.

NB. colorings.
input =: >;: 'pgpygr rprrgy ppryyg rrygpy' 

NB. all 24 rotations of a cube.
gens =: (i.6),'NSKTFL'i.'FSLTKN',:'SKTNFL'
rots =: ([:~.[:/:~[:,/{"1/)^:_~ gens
   
NB. all 24 rotations of each 4 cube.
posn =: rots {"_ 1 input

NB. stack a list of 4 lists of 24 cubes to a tower of 4 cubes in each 
+24^4 way,
NB. order of 4 cubes doesn't matter.
stack =: [:>[:>[:,[:{[:<"1<"1
towers =: stack posn

NB. check if a cube is good or bad -- rank 2, gives bool atom.
good =: [: (*./"1) 4 ({."1) 4 = [: (#@:~."1) |:"2
solns =: (good # ]) towers

NB. each solution gives 8 trivially equivalent solutions because you c
+an 
NB. rotate all cubes simultanously around the vertical axis or twice a
+round
NB. the horizontal axis.  (this may be less than 8 for degenerate cube
+s.)
NB. rota y gives all 8 rotations of a tower y
rott =: (2{gens)&({"_ 1)
rotu =: ({~1{gens)&({"_ 1)
rota =: ([: ,/ [: rott^:(<4) rotu^:(<2))"2
NB. verify that each rotation of each solution is indeed among the sol
+utions.
assert =: 3 :'assert. y'
assert *./ solns e.~ ,/ rota solns
NB. now normalize the solutions by changing a tower to the 
NB. lexicographically first out of these eight variants,
NB. and weed out repetitions.
norml =: ([:{.[:/:~rota)"2
solns1 =: ~./:~ norml solns

NB. output the solutions. there's just one if you use the sample
NB. cube colorings in the post.
echo solns1
echo #solns1

NB. end
[download]