Polya Counting

How many different ways are there to colour the faces of a dodecahedron using red, green and blue? No, the answer's not 312 - two colourings are equivalent if you can turn one into the other by rotating the dodecahedron - and what we want to count is the number of equivalence classes. If you spend a little time thinking about it, you'll quickly see that there are an awful lot, and that enumerating them all is no easy task.

However, there is a trick we can use - Polya counting - and that is the subject of this article. To explain how it works, let's start with a simpler example. How many ways are there to colour the faces of a cube with two colours, up to rotation? With a little thought, you will see that there are 10, as shown below:

But we don't want to have to think about it, we want to be able to just turn a handle and get the answer. So how do we do it? Well, first of all, we need to find the rotation group of the cube. We'll number the faces like this:

(So I'm thinking of 1 as the top, 6 as the bottom, and 2,3,4,5 as the sides.)

We need to find generators for the group. We'll try a face rotation, a vertex rotation, and an edge rotation: 

import PermutationGroups
import SchreierSims
import MPoly
import QQ

cubeF = fromCycles [[1],[2,3,4,5],[6]]
cubeV = fromCycles [[1,2,3],[4,5,6]]
cubeE = fromCycles [[1,3],[2,4],[5,6]]

cubeGp = schreierSimsTransversals 6 [cubeF, cubeV, cubeE]

> orderSS cubeGp
24

This is the right answer, so we have the full group. Let's have a look at the elements. 

> mapM_ print (eltsSS cubeGp)
[[1,4,6,2]]
[[1,3,2],[4,6,5]]
[[1,5,2],[3,4,6]]
[[1,2],[3,5],[4,6]]
[]
[[1,3,6,5]]
[[1,5,6,3]]
[[1,6],[3,5]]
[[1,2,3],[4,5,6]]
[[2,3,4,5]]
[[1,6],[2,3],[4,5]]
[[1,4,5],[2,3,6]]
[[1,3],[2,4],[5,6]]
[[2,4],[3,5]]
[[1,6],[2,4]]
[[1,5],[2,4],[3,6]]
[[2,5,4,3]]
[[1,2,5],[3,6,4]]
[[1,4,3],[2,5,6]]
[[1,6],[2,5],[3,4]]
[[1,3,4],[2,6,5]]
[[1,4],[2,6],[3,5]]
[[1,2,6,4]]
[[1,5,4],[2,6,3]]

We see that there are five different types of elements:

Now, Polya counting relies on a clever way of writing this information as a multivariate polynomial, called the cycle index. It's easiest to define by example:

You get the idea. The cycle index of an element comes from taking the cycle type and replacing each n-cycle by the variable xn.

Next, the cycle index of the whole group is the sum of the cycle indices of all the elements, divided by the order of the group. In this case
1/24 (6 x12 x41 + 3 x12 x22 + 8 x32 + 6 x23 + x16) .

In Haskell: 

cycleIndexElt g = monomial (cycleType g)
cycleIndexGp gp = sum [cycleIndexElt g | g <- eltsSS gp] / fromInteger (orderSS gp)

> cycleIndexElt cubeF
x1^2x4
> cycleIndexGp cubeGp
1/24x1^6+1/8x1^2x2^2+1/4x2^3+1/4x1^2x4+1/3x3^2

Now comes the good bit. We have two colours, blue and red, which we now equate with the variables x and y. We make the following substitutions:

When we make these substitutions into the cycle index, we get what I'll call the counting polynomial. In Haskell: 

countingPoly gp vs = substMP (cycleIndexGp gp) [sum [v^i | v <- vs] | i <- [0..]]

> countingPoly cubeGp [x,y]
x^6+x^5y+2x^4y^2+2x^3y^3+2x^2y^4+xy^5+y^6

Now for the amazing bit. This polynomial is telling us something about the colourings of the cube. Each monomial corresponds to a colouring of the cube. The x6 monomial tells us that there is one colouring with 6 blue faces. The 2 x4 y2 tells us that there are two colourings with 4 blue and 2 red faces. And so on. Here are the cube colourings again, with the polynomial written underneath.

So to answer our question - to count the total number of colourings - we just need to sum the coefficients of the counting polynomial.

What's going on? Why does this work? To understand that, we're going to need to "remember" some of the structure that we "forgot" when we moved from coloured cubes and rotation groups to multivariate polynomials.

First of all, let's think about the cycle index again. Consider the face rotation [[1],[2,3,4,5],[6]]. Its cycle index is x12x4, meaning that it is composed of two 1-cycles and a 4-cycle. Let's look at the 4-cycle first. What happens when we make the substitution x4 = x4 + y4 ? Remember that x means red, y means blue. x4 = x4 + y4 is saying that if you have a 4-cycle, then you can colour it either as four reds or as four blues. Anything else - a mixture of reds and blues - and the colouring won't be preserved by this group element. The two substitutions x1 = x + y are saying (trivially) that each of the 1-cycles can be coloured either red or blue. Then, the fact that we multiply these factors together is saying that these choices can be made independently of one another.

When we make the substitutions and multiply out, we get
x12 x4 = (x+y)2 (x4 + y4)
= x6 + 2 x5 y + x4 y2 + x2 y4 + 2 x y5 + y6
The terms in the result are counting the colourings of the cube that stay the same under the action of this group element.

So let's consider what we've done. Our group elements act naturally on the space of cube faces, numbered 1 to 6 - call it X. However, we've now switched to considering it's action on a different space, the space of all 26 possible colourings of the cube - call it Y. The counting polynomial of a group element - the polynomial you get when you make the substitutions into the cycle index of an element - is counting the number of points in this new space Y which are fixed by the group element.

Now, the cycle index of the group was defined as the "average" of the cycle indices of the elements. So when we make the substitution, the counting polynomial of the group will be the average of the counting polynomials of the elements - which count the number of points fixed by the elements.

But now Burnside's lemma states that when a group acts on a space, then the number of orbits is equal to the "average" number of fixed points of its elements.
|X/G| = sum [fix g | g <- G] / |G|
Which means we're done. Burnside's lemma says that the number of orbits in the space Y of cube colourings is equal to the average number of fixed points - which is what the cycle index counts. But the number of orbits in Y is precisely the thing we were looking for - the number of colourings modulo rotational equivalence.

Whew! That took rather longer than I intended - but we made it. (Well, it probably took me longer to write it than it took you to read it.)

Now, in case you'd forgotten, the reason we did all that is because we want to find out how many ways we can colour a dodecahedron with three colours, up to rotation. Well we should be able to whizz through it now. Let's start by finding the rotation group of a dodecahedron. We'll number the faces of the dodecahedron as follows:

(Oops - this is the dodecahedron before it's been stuck together. On the left are the top six pentagons, on the right the bottom six, both sets seen from above. We make the dodecahedron by sticking the top and bottom together. When we do that, the numbers on opposite faces will add up to 13 - check that this is the way you are visualising it.)

Now, as before, let's write down some generators: 

dodecF = fromCycles [[1],[2,3,4,5,6],[11,10,9,8,7],[12]]       -- rotation about a face
dodecV = fromCycles [[1,2,3],[4,6,8],[9,7,5],[12,11,10]]       -- rotation about a vertex
dodecE = fromCycles [[1,2],[3,6],[4,9],[5,8],[7,10],[11,12]]   -- rotation about an edge

dodecGp = schreierSimsTransversals 12 [dodecF, dodecV, dodecE]

> orderSS dodecGp
60

A quick internet search confirms that this is the order of the full rotation group - eg here. Let's carry on: 

> cycleIndexGp dodecGp
1/60x1^12+1/4x2^6+1/3x3^4+2/5x1^2x5^2

Now, this time round we have three colours, which we equate with the variables x, y and z. So we need to make the following substitutions: 

> countingPoly dodecGp [x,y,z]
x^12+x^11y+3x^10y^2+5x^9y^3+12x^8y^4+14x^7y^5+24x^6y^6+14x^5y^7+12x^4y^8+5x^3y^9
+3x^2y^10+xy^11+y^12+x^11z+3x^10yz+11x^9y^2z+33x^8y^3z+66x^7y^4z+94x^6y^5z+94x^5
y^6z+66x^4y^7z+33x^3y^8z+11x^2y^9z+3xy^10z+y^11z+3x^10z^2+11x^9yz^2+57x^8y^2z^2+
132x^7y^3z^2+246x^6y^4z^2+278x^5y^5z^2+246x^4y^6z^2+132x^3y^7z^2+57x^2y^8z^2+11x
y^9z^2+3y^10z^2+5x^9z^3+33x^8yz^3+132x^7y^2z^3+312x^6y^3z^3+462x^5y^4z^3+462x^4y
^5z^3+312x^3y^6z^3+132x^2y^7z^3+33xy^8z^3+5y^9z^3+12x^8z^4+66x^7yz^4+246x^6y^2z^
4+462x^5y^3z^4+600x^4y^4z^4+462x^3y^5z^4+246x^2y^6z^4+66xy^7z^4+12y^8z^4+14x^7z^
5+94x^6yz^5+278x^5y^2z^5+462x^4y^3z^5+462x^3y^4z^5+278x^2y^5z^5+94xy^6z^5+14y^7z
^5+24x^6z^6+94x^5yz^6+246x^4y^2z^6+312x^3y^3z^6+246x^2y^4z^6+94xy^5z^6+24y^6z^6+
14x^5z^7+66x^4yz^7+132x^3y^2z^7+132x^2y^3z^7+66xy^4z^7+14y^5z^7+12x^4z^8+33x^3yz
^8+57x^2y^2z^8+33xy^3z^8+12y^4z^8+5x^3z^9+11x^2yz^9+11xy^2z^9+5y^3z^9+3x^2z^10+3
xyz^10+3y^2z^10+xz^11+yz^11+z^12

Hmm, looks like there are quite a few. Let's write a little more code, to sum up the coefficients. 

polyaCount gp vs = sum [c | (c,as) <- termsMP (countingPoly gp vs)]

> polyaCount dodecGp [x,y,z]
9099

So there you are - quite a few - good thing we didn't try to count them by hand.

Next time, I'll be looking at how to use Polya counting to count the number of non-isomorphic graphs on a given number of vertices. Incidentally, all the code we wrote here is in PolyaCounting.hs .

Copyright (c) David Amos, 2006
polyomino (at) f2s (dot) com