1.5 Enumerating graphs with Pólya’s theorem and GMP

by Marko Riedel

1.5.1 Introduction

The object of this recipe is to enumerate non-isomorphic graphs on n vertices using Pólya’s theorem and GMP (the GNU multiple precision arithmetic library). We will be concerned with the algorithmics rather than the mathematics. You are encouraged to learn more about the latter, e.g. from Frank Harary’s books on the subject.

Fig. 1 shows all possible graphs on three vertices. There are eight of these; for a graph on n vertices we have 2^n(n-1)∕2 such graphs. We want a different statistic, namely one that takes isomorphisms into account. Fig. 2 shows that there are really only four different graphs on three vertices. This is because e.g. there exists a vertex permutation that maps the three graphs with one edge to one another. We say that these three graphs constitute an equivalence class.

Let us now reformulate the problem so that Pólya’s theorem applies. Given n, we want to know how many non-isomorphic graphs, i.e. equivalence classes there are with k edges. We can think of the edges as slots that are either filled (edge present) or empty (edge absent). An empty slot has size zero and a filled one, size one. The number of edges is the sum of the sizes of the slots. The slots (edge locations) are being permuted by any of the n! vertex permutations. There is one edge permutation for each vertex permutation. We will obtain the values for all k at once, namely in a generating function, i.e. in a function where the term qz^k denotes the presence of q non-isomorphic graphs with k edges.

The generic definition of our problem is this: given a set of slots that are being permuted by a permutation group (more on this later) and a set of objects of different sizes that go into the slots (one object per slot), count the number of equivalence classes of filled slot configurations with respect to the permutation group whose total size (i.e. sum of the sizes of the objects placed in the slots) is k. In our case the group is the pair group obtained from the symmetric group acting on n vertices and the objects that go into the slots are present/absent edges.

Example. Recall Fig. 1 and label the vertices 1, 2, 3. The table shows the vertex and the corresponding edge permutations (the three edges are 12, 13 and 23; the permutations are in table notation, e.g. 1 → 2, 2 → 3 and 3 → 1 is written as 231):

The eight graphs correspond to the three slots (edges) being filled or empty.

The edge permutations permute the columns of the table, e.g. 12 23 13 takes line two to line three, so these are in the same equivalence class. Permuting edges does not change the total size of the configuration (number of edges), so there is no permutation that takes line two to line four.

Another common application of Pólya’s theorem is to count, say, the number of different necklaces, where a necklace is made up of n gemstones and there are m different types of gemstones (colors) available. We can choose these to have the same size or different sizes. In terms of the above discussion, the slots are the positions on the necklace, the objects that go into the slots are the gemstones and the permutations rotate or flip the necklace or both.

1.5.2 Permutations

We will be using only a few basic facts about permutations and permutation groups. The elements of the latter are permutations and the multiplication operation is permutation composition. There must not be any permutations whose product is not in the group. E.g. {123,231} is not a group because 231 ∘ 231 = 312, which is not in the set. We are concerned with two groups: the symmetric group S_n on n vertices, which includes all n! permutations and the corresponding pair group S_n⁽²⁾ on n(n - 1)∕2 edges, whose elements are obtained from the action of the vertex permutation on the edges. E.g. consider the vertex permutation 23451 and the edges 12 13 14 15 23 24 25 34 35 45. The corresponding edge permutation ε is 23 24 25 12 34 35 13 45 14 15.

We now switch to a more useful way of writing permutations, i.e. the so-called disjoint cycle decomposition. This consists in writing the permutation as a product of its cycles. E.g. the permutation 164325 (table notation) would be written as (1) (265) (34) since e.g. it takes 2 to 6, 6 to 5 and 5 back to 2, and the edge permutation ε as (12 23 34 45 15) (13 24 35 14 25).

The disjoint cycle decomposition tells us about the shape of the permutation. This shape information is captured by a product notation, where the variable a_k stands for a cycle of length k and j_k for the number of cycles with this length, giving the product ∏ _k=1^pa_k^j_k. E.g. the product for (1) (265) (34) would be a ₁a₂a₃, whith j₁ = j₂ = j₃ = 1 and j₄ = j₅ = j₆ = 0 and for ε, a₅², with j₁ = j₂ = j₃ = j₄ = 0, j₅ = 2 and j₆ = j₇ = j₈ = j₉ = j₁₀ = 0.

1.5.3 The Cycle Index

The cycle index Z(S) of a permutation group S is the average of ∏ _k=1^pa_k^j_k over all permutations in the group. E.g. the cycle indices of the symmetric group S₃ and the pair group S₃⁽²⁾ are obtained from the following table.

This special case yields

1.5.4 Pólya’s theorem

This is what the theorem says: suppose you have some set of objects of different sizes and distribute these objects into n slots, where a permutation group acts on the slots to create equivalence classes of filled slot configurations. The generating function of these equivalence classes is obtained by replacing a_k by f(z^k) in the cycle index, where f is the generating function of the objects.

The permutation group is S_n⁽²⁾ (pair group acting on the vertices) when we enumerate graphs and the generating function of the objects is 1 + z, indicating whether an edge is present (size one, z) or not (1). Substitute a_k = 1 + z^k into Z(S₃⁽²⁾) to obtain

which says that there is one graph with three edges, one with two, one with one edge and one with no edges, as is also evident from Fig.2.

1.5.5 Cycle index of the pair group

We need the cyle index of S_n⁽²⁾ in order to enumerate graphs. There is one permutation for every permutation from S_n in S_n⁽²⁾. Fortunately we needn’t iterate over all n! permutations in S_n, which would make the computation infeasible for all but small n. Recall that the shape of the permutation determines the product that it contributes to the cycle index. This means that e.g. (1 2 3) (4) (5 6), (1) (2 3) (4 5 6) and (1 2) (3 4 5) (6) all contribute the same term to S₆, namely a₁a₂a₃ (three copies of it).

Hence we can enumerate the permutations in S_n by iterating over all possible shapes and computing how often each shape occurs. These shapes are precisely the partitions of n into some number of summands. E.g. the cycle index for S₃ contains one product term for each of the possible partitions 1 + 1 + 1, 1 + 2 and 3. We need the coefficient associated to each of these three partitions, i.e. the number of times this shape occurs.

The disjoint cycle decomposition ∏ _k=1ⁿa_k^j_k occurs

times. (Partition the n elements into subsets, one for each cycle. A subset of size k generates k!∕k cycles. The same size j_k set of k-cycles may be chosen in j_k! different ways.)

For S₃ this yields the following results:

• 1+1+1: choose three subsets of size one in 3!∕1!∕1!∕1! = 6 ways, every such subset generates 1!∕1 = 1 cycles but all 3! choices for the three subsets are the same, giving 1.
• 1 + 2: choose a subset of size two and one of size one in 3!∕1!∕2! = 3 ways (there are three ways to choose the singleton): the size two subset generates 2!∕2 = 1 cycle and the size one subset, 1!∕1 = 1 cycle; there no duplicates in the choice of size 2 or size 1 subsets, giving 3.
• 3: choose a subset of size three in 3!∕3! = 1 ways, this subset generates 3!∕3 = 2 cycles and again there are no duplicates in the choice of the size 3 subset, giving 2.

This computation confirms the value for S₃ that we computed earlier.

How to obtain the possible partitions? This is simple, since there is a folklore algorithm for the problem. Instead of computing the set of partitions P_n of n, we start with P_n,k, the set of partitions of n into k summands. Now every partition in P_n,k either contains one, and can be obtained by including one in a partition from P(n- 1,k - 1), or it does not contain one and n-k >= k, and it can be obtained from a partition from P(n-k,k) by incrementing its elements. We can compute P_n,k recursively, the base case being P_1,1 = {[1]}. P_n is the union of the P_n,k for 1 ≤ k ≤ n.

At this point we have all possible shapes and we know how often they occur. The last step is to compute the shape of the edge permutation from the shape of the vertex permutation. An edge connects two vertices. There are three possiblities: these vertices lie on the same cycle of the vertex permutation, they lie on different cycles of the same length or they lie on different cycles of different lengths. These three possibilities are shown in Fig. 3. We must compute the length of the edge cycle (i.e. how long it takes for the vertex pair to return to the start edge when it moves in parallel along the vertex cycle) in each case and how often it occurs. There is an additional distinction to be made when the vertices lie on the same cycle: a cycle of even length contains vertex pairs that return to the start edge after having covered exactly half, rather than the whole cycle. Furthermore, the number of steps for the return to the start position for two vertices on cycles of different lengths is given by the least common multiple of those lengths. Note that we’ll overcount cycles by a factor that is equal to the length of the edge cycle. We illustrate with a sample computation.

Say the vertex permutation has the shape a₃² a₆ a₇.

• Same cycle.
  • a₃²: three vertex pairs, length three, twice: a₃²
  • a₆: fifteen vertex pairs, three short cycles: a₆²a₃
  • a₇: twenty-one vertex pairs, length seven, a₇³
• Different cycles, same length.
  • a₃²: nine vertex pairs, length three: a₃³
• Different cycles, different length.
  • a₃²a₆: eighteen vertex pairs, length six, twice, a₆⁶
  • a₃²a₇: twenty-one vertex pairs, length twenty-one, twice, a₂₁²
  • a₆a₇: forty-two vertex pairs, length forty-two: a₄₂.

This shows that the corresponding edge permutation has the shape a₃⁶ a₆⁸ a₇³ a₂₁² a₄₂, and indeed 18 + 48 + 21 + 42 + 42 = 171 = 1∕2 19 (19 - 1).

We now have all the ingredients to implement the graph enumeration algorithm. There will be four steps:

• Compute the partitions of n.
• Compute the vertex permutation shapes from the partitions and the coefficient.
• Compute the edge permutation shape, substitute 1+z and add the resulting polynomial to the current total.
• Divide by n! to obtain the generating function.

Note that we average the cycle index in the last step in order to be able to work with integers during the whole computation. Speaking of integers, our implementation will use the GMPInt wrapper class for GMP integers.

1.5.6 Concerning memory allocation

The computation that we undertake produces many intermediate results. These are: GMPInt values to represent coefficients of the products whose sum forms the cycle index; NSString objects that document the computation if the user activates debugging and NSMutableArray objects that hold partitions and product terms as well as polynomials in one variable that result from the substitution of 1 + z into the cycle index. All of these can use up a lot of memory very quickly, which is why we absolutely must implement garbage collection. This is done in the loop that iterates over the partitions and hence over the products that occur in the cycle index. We check whether the number of allocated critical objects has doubled, and empty the autorelease pool if it has.

1.5.7 Implementation I: auxiliary functions

There is a preliminary remark concerning our terminology. We use the term multiset to refer to partitions and products. This is because a partition like [1,1,1,2,3] or a₁³ a₂ a₃ is a multiset containing three copies of the value “1” and one copy each of the values “2” and “3.”

Our list of auxiliary functions starts with the basic GCD algorithm. The GCD g of p and q divides p = mq + r and q, so it also divides r < q. Hence replace (p,q) by (q,r) and repeat until r = 0. We need the GCD to compute the LCM of the lengths of two cycles of a vertex permutation. The LCM gives the length of the edge cycle of those edges that have one vertex on each of the two vertex cycles.

The function SumIt sums a multiset, where each element contributes according to its multiplicity and the first element has value one. It iterates over the elements in a loop. This is very useful for consistency checks. The cycle lengths of a vertex permutation must sum to n, and the cycle lengths of the corresponding edge permutation to 1∕2n(n - 1).

The next two functions are very important; they compute P_n,k and P_n, respectively. Start with the former. The base case comes first: the answer is {[1]} when n = 1. The data structure that we use is an array of arrays of integer NSNumber objects. The number is one summand from the partition, the inner array a single partition and the outer array a collection of partitions.

The two recursions are next. Either the partition contains one or it doesn’t. There is a one-to-one correspondence with a permutation of P_n-1,k-1 in the former case, so obtain those and collect them, adding a summand with value one at the front of each summand array.

If the partition does not contain one and n - k is large enough for k summands, then we have a one-to-one correspondence with a permutation from P_n-k,k. We obtain these permutations and increment each summand of each partition, collecting the new partitions in the array result. These are the only possible cases and we are done.

Obtaining all partitions of n is now easy: compute the union of P_n,k for 1 ≤ k ≤ n. This is done with a loop.

We will need to turn partitions into multisets, i.e. vertex permutation shapes that represent product terms of the form ∏ _k=1ⁿa_k^j_k. This product correponds to a partition that contains j _k instances of the value k. E.g. we want to turn [1,1,1,2,2,3] into a₁³ a₂² a₃, which is represented as the array [3,2,1]. First initialize the array with zero, i.e. in our example we would have [0,0,0] and the capacity of the array would be three.

Then iterate over the partion and record a summand with value value (recall that these values are wrapped in NSNumber objects) at position value-1 of the multiset array. This is done by reading the old value, incrementing it, and writing it back into the array. This would yield the sequence [1,0,0], [2,0,0], [3,0,0], [3,1,0], [3,2,0], and [3,2,1] in the example.

Our program provides an extensive trace/debug facility that makes it possible for the user to follow the computation very closely. Hence it is imperative that we be able to convert from the internal representation of a multiset into a string representation, i.e. we must be able to turn [3,2,1] into the string a₁³ a₂² a₃. This is done by iterating over the multiset, reading the entry at that position, and storing a string representation of the term in the array factors. (Empty fields are skipped.)

The string representation of such a term has the form a_k if if j_k is one, and the form a_k^j_k otherwise.

We join all terms that have been collected by concatenating them. (They are separated by a single space in the result string.)

Recall that the number of times a shape occurs in the symmetric group is given by the formula

The next procedure computes the denominator of this fraction. Set up a loop to iterate over the multiset and process nonzero entries only.

The actual computation is to multiply the current denominator by k^j_k and j _k!.