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 2n(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.


PIC

Figure 1: All eight graphs on three vertices.



PIC

Figure 2: All non-isomorphic graphs on three vertices.


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 qzk 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):



1 2 312 13 23


1 3 213 12 23


2 1 312 23 13


2 3 123 12 13


3 1 213 23 12


3 2 123 13 12


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




12 13 23






emptyemptyempty



emptyempty filled



empty filled empty



empty filled filled



filled emptyempty



filled empty filled



filled filled empty



filled filled filled



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 Sn on n vertices, which includes all n! permutations and the corresponding pair group Sn(2) 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 ak stands for a cycle of length k and jk for the number of cycles with this length, giving the product k=1pakjk. E.g. the product for (1) (265) (34) would be a 1a2a3, whith j1 = j2 = j3 = 1 and j4 = j5 = j6 = 0 and for ε, a52, with j1 = j2 = j3 = j4 = 0, j5 = 2 and j6 = j7 = j8 = j9 = j10 = 0.

1.5.3 The Cycle Index

The cycle index Z(S) of a permutation group S is the average of k=1pakjk over all permutations in the group. E.g. the cycle indices of the symmetric group S3 and the pair group S3(2) are obtained from the following table.





1 2 3 a13 12 13 23 a13




1 3 2a1a213 12 23a1a2




2 1 3a1a212 23 13a1a2




2 3 1 a3 23 12 13 a3




3 1 2 a3 13 23 12 a3




3 2 1a1a223 13 12a1a2




This special case yields

                1     1      1
Z(S3) = Z(S(23)) =-a31 +-a1a2 +- a3.
                6     2      3

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 ak by f(zk) in the cycle index, where f is the generating function of the objects.

The permutation group is Sn(2) (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 ak = 1 + zk into Z(S3(2)) to obtain

1(1+z)3+ 1(1+z)(1+z2)+  1(1+z3) = 1z3+ 1z2+ 1z+ 1+ 1z3+ 1z2+ 1z +1 + 1z3+ 1= z3+z2 +z +1,
6        2              3         6    2    2   6  2    2    2   2   3    3

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 Sn(2) in order to enumerate graphs. There is one permutation for every permutation from Sn in Sn(2). Fortunately we needn’t iterate over all n! permutations in Sn, 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 S6, namely a1a2a3 (three copies of it).

Hence we can enumerate the permutations in Sn 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 S3 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=1nakjk occurs

          ∏n (   )jk∏n
∏---n!----     k!      -1-=  ∏n--n!----
  nk=1(k!)jkk=1  k    k=1jk!    k=1 kjkjk!

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

For S3 this yields the following results:

This computation confirms the value for S3 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 Pn of n, we start with Pn,k, the set of partitions of n into k summands. Now every partition in Pn,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 Pn,k recursively, the base case being P1,1 = {[1]}. Pn is the union of the Pn,k for 1 k n.


PIC

Figure 3: Vertex permutation with the shape a43 a8 and the three possible types of edges.


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 a32 a6 a7.

This shows that the corresponding edge permutation has the shape a36 a68 a73 a212 a42, and indeed 18 + 48 + 21 + 42 + 42 = 171 = 12 19 (19 - 1).

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

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 a13 a2 a3 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.

 
int gcd(int p, int q) 
{ 
    int r; 
 
    if(p<q){ 
        int s; 
 
        s = p; 
        p = q; 
        q = s; 
    } 
 
    while(r=p%q){ 
        p = q; 
        q = r; 
    } 
 
    return q; 
}

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 12n(n - 1).

 
int SumIt(NSArray *data) 
{ 
    int sum, pos, max = [data count]; 
 
    for(sum=0, pos=0; pos<max; pos++){ 
        sum += (pos+1)* 
            [[data objectAtIndex:pos] intValue]; 
    } 
 
    return sum; 
}

The next two functions are very important; they compute Pn,k and Pn, 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.

 
NSMutableArray *Partitionsk(int n, int k) 
{ 
    NSMutableArray *part, 
        *result = [NSMutableArray arrayWithCapacity:1]; 
    NSEnumerator *en; 
 
    if(n==1){ 
        part = 
            [NSMutableArray 
                arrayWithObject: 
                     [NSNumber numberWithInt:1]]; 
        [result addObject:part]; 
        return result; 
    }

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 Pn-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(k>1){ 
        en = [Partitionsk(n-1, k-1) objectEnumerator]; 
        while((part = [en nextObject])!=nil){ 
            [part insertObject:[NSNumber numberWithInt:1] 
                   atIndex:0]; 
            [result addObject:part]; 
        } 
    }

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 Pn-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.

 
    if(n-k>=k){ 
        en = [Partitionsk(n-k, k) objectEnumerator]; 
        while((part = [en nextObject])!=nil){ 
            NSMutableArray *npart = 
                [NSMutableArray arrayWithCapacity:k]; 
            int pos; 
            for(pos=0; pos<k; pos++){ 
                int val = [[part objectAtIndex:pos] intValue]; 
                [npart 
                     addObject: 
                         [NSNumber numberWithInt:val+1]]; 
            } 
            [result addObject:npart]; 
        } 
    } 
 
    return result; 
}

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

 
NSMutableArray *Partitions(int n) 
{ 
    int k; 
    NSMutableArray *result = 
        [NSMutableArray arrayWithCapacity:1]; 
 
    for(k=1; k<=n; k++){ 
        NSMutableArray *partsk = Partitionsk(n, k); 
        [result addObjectsFromArray:partsk]; 
    } 
 
    return result; 
}

We will need to turn partitions into multisets, i.e. vertex permutation shapes that represent product terms of the form k=1nakjk. 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 a13 a22 a3, 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.

 
NSMutableArray *PartitionToMSet(NSMutableArray *part) 
{ 
    int pos, max = [[part lastObject] intValue]; 
    NSMutableArray *result = 
        [NSMutableArray arrayWithCapacity:max]; 
 
    for(pos=0; pos<max; pos++){ 
        [result addObject:[NSNumber numberWithInt:0]]; 
    }

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.

 
    NSEnumerator *en = [part objectEnumerator]; 
    NSNumber *val; 
    while((val = [en nextObject])!=nil){ 
        int value = [val intValue], 
            count = [[result objectAtIndex:value-1] intValue]; 
        [result replaceObjectAtIndex:value-1 
                withObject:[NSNumber numberWithInt:count+1]]; 
    } 
 
    return result; 
}

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 a13 a22 a3. 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.)

 
NSString *MSetToProduct(NSMutableArray *mset) 
{ 
    NSNumber *cval; 
    NSMutableArray *factors = 
        [NSMutableArray arrayWithCapacity:1]; 
    int pos, max = [mset count]; 
 
    for(pos=0; pos<max; pos++){ 
        int count = [[mset objectAtIndex:pos] intValue];

The string representation of such a term has the form ak if if jk is one, and the form akjk otherwise.

 
        if(count){ 
            NSString *item = 
                [NSString stringWithFormat:@”a_%d”, 
                           pos+1]; 
            if(count>1){ 
                item = [item stringByAppendingFormat:@”ˆ%d”, 
                              count]; 
            } 
            [factors addObject:item];

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

 
        } 
    } 
 
    return [factors componentsJoinedByString:@”_”]; 
}

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

    n!
∏n---kjkj!.
  k=1    k

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

 
GMPInt *MSetToCoeffDenom(NSMutableArray *mset) 
{ 
    GMPInt *one = [GMPInt oneObj], 
        *res = one, *k = one; 
    int pos, max = [mset count]; 
 
    for(pos=0; pos<max; pos++){ 
        int count = [[mset objectAtIndex:pos] intValue]; 
        if(count){

The actual computation is to multiply the current denominator by kjk and j k!.

 
            GMPInt 
                *m1 = [k powUI:count], 
                *m2 = [GMPInt factorialUI:count]; 
 
            res = [res mul:m1]; 
            res = [res mul:m2]; 
        }

We increment k at the end of the loop body and return the denominator once the loop exits.

 
        k = [k add:one]; 
    } 
 
    return res; 
}

The next function is really more of a macro that we’ll use to turn vertex permutation shapes into edge permutation shapes, which is why it’s marked inline. It records an increase of a product term ak by some amount.

 
inline void AddToMSet(NSMutableArray *mset, int len, int count) 
{ 
    int val = [[mset objectAtIndex:len-1] intValue]; 
    [mset replaceObjectAtIndex:len-1 
          withObject:[NSNumber numberWithInt:val+count]]; 
}

The next function is of critical importance: it produces the terms of the cycle index of the pair group. The maximum possible length of an edge cycle obtained from a vertex permutation is l(l - 1), where l is the length of the longest cycle of the vertex permutation. This occurs if there are cycles of length l - 1 and applies to edges with one vertex on the l - 1-cycle and the other on the l-cycle. The maximum length is one for the identity permutation. The multiset of edge cycles is initialized to be empty.

 
NSMutableArray *VertexMSetToEdgeMSet(NSMutableArray *vmset) 
{ 
    int vmxlen = [vmset count], vlen1, vlen2, 
        upper = 0, 
        emxlen = (vmxlen==1 ? 1 : vmxlen*(vmxlen-1)), elen; 
    int inst1, inst2, instcount; 
    NSMutableArray *emset = 
        [NSMutableArray arrayWithCapacity:emxlen]; 
 
    for(elen=1; elen<=emxlen; elen++){ 
        [emset addObject:[NSNumber numberWithInt:0]]; 
    }

The first case occurs when there are two vertices on the same cycle. If the cycle length is odd, then the edge cycle has the same length and the number of edges on such cycles is the number of different vertex pairs on the cycle. We scale (divide) by the length of the cycle because each cycle is counted once for each edge on it. We must also scale the term by its multiplicity. (This corresponds to choosing one of jk cycles.) We also record the new maximum edge cycle length if it has changed. The latter two operations, i.e. compensating for overcount and recording the new maximum edge cycle length are common to all cases and will not be discussed further.

 
    // two vertices on the same cycle 
    for(vlen1=2; vlen1<=vmxlen; vlen1++){ 
        inst1 = [[vmset objectAtIndex:vlen1-1] intValue]; 
        if(inst1){ 
            if(vlen1%2){ 
                elen = vlen1; 
                instcount = vlen1*(vlen1-1)/2; 
                instcount *= inst1; 
 
                instcount /= elen; 
                AddToMSet(emset, elen, instcount); 
                if(elen>upper){ 
                     upper = elen; 
                } 
            }

The next case corresponds to two vertices being located on a cycle of even length. Edges whose vertices are directly accross from one another need only move through half the cycle to return to the start position. This case is handled in the first half of the clause. The second half corresponds to vertex pairs that are not directly accross from one another and need to move through the entire cycle.

 
            else{ 
                elen = vlen1/2; 
                instcount = vlen1/2; 
                instcount *= inst1; 
 
                instcount /= elen; 
                AddToMSet(emset, elen, instcount); 
                if(elen>upper){ 
                     upper = elen; 
                } 
 
                elen = vlen1; 
                instcount = vlen1*(vlen1-1)/2-vlen1/2; 
                instcount *= inst1; 
 
                instcount /= elen; 
                AddToMSet(emset, elen, instcount); 
                if(elen>upper){ 
                     upper = elen; 
                } 
 
            } 
        } 
    }

The second case occurs when two vertices lie on different cycles that have the same length. We can choose these two cycles in 12jk(jk - 1) ways. The two vertices move in parallel along the two cycles and return to the start point only after having completed the whole vertex cycle.

 
    // two vertices on different cycles of the same length 
    for(vlen1=1; vlen1<=vmxlen; vlen1++){ 
        inst1 = [[vmset objectAtIndex:vlen1-1] intValue]; 
        if(inst1>1){ 
            elen = vlen1; 
            instcount = vlen1*vlen1; 
            instcount *= inst1*(inst1-1)/2; 
 
            instcount /= elen; 
            AddToMSet(emset, elen, instcount); 
            if(elen>upper){ 
                upper = elen; 
            } 
        } 
    }

The last case occurs when the two vertices lie on different cycles of different lengths. We check all possible combinations. We must process all pairs (jk1,jk2) where jk1 and jk2 are not zero.

 
    // two vertices on different cycles of different lengths 
    for(vlen1=1; vlen1<=vmxlen; vlen1++){ 
        for(vlen2=vlen1+1; vlen2<=vmxlen; vlen2++){ 
            inst1 = [[vmset objectAtIndex:vlen1-1] intValue]; 
            inst2 = [[vmset objectAtIndex:vlen2-1] intValue]; 
 
            if(inst1 && inst2){

The number of steps is the LCM of the lengths of the two cycles. E.g. for a cycle of length 2 and one of length 3, we visit 11, 22, 13, 21, 12, 23 and 11, a total of six steps. There are k1k2 possible edges for each of jk1jk2 pairs of vertex cycles.

 
                elen = vlen1*vlen2/gcd(vlen1, vlen2); 
                instcount = vlen1*vlen2; 
                instcount *= inst1*inst2; 
 
                instcount /= elen; 
                AddToMSet(emset, elen, instcount); 
                if(elen>upper){ 
                     upper = elen; 
                } 
            } 
        } 
    }

The edge permutation shape is now ready. We remove any zeros that may remain at the upper end.

 
    if(upper<emxlen){ 
        [emset 
            removeObjectsInRange: 
                NSMakeRange(upper, emxlen-upper)]; 
    } 
 
    return emset; 
}

The remaining functions implement basic arithmetic for polynomials in one variable. We use algorithms that suffice for our needs. E.g. we do not implement polynomial multiplication by the FFT.

1.5.8 Implementation II: Polynomials

Polynomials are represented by arrays of GMPInt coefficients. Start with a very basic polynomial, namely the term f(zk), i.e. 1 + zk. We set all coefficients except for the constant term and zk to zero.

 
NSMutableArray *PolyK(int k) 
{ 
    NSMutableArray *p 
        = [NSMutableArray arrayWithCapacity:k+1]; 
    int l; 
    GMPInt *zero = [GMPInt zeroObj], *one = [GMPInt oneObj]; 
 
    [p addObject:one]; 
    for(l=1; l<k; l++){ 
        [p addObject:zero]; 
    } 
    [p addObject:one]; 
 
    return p; 
}

We add polynomials by iterating over their coefficients and computing the sum of the coefficients of like powers.

 
NSMutableArray *PolySum(NSMutableArray *p1, NSMutableArray *p2) 
{ 
    int c1 = [p1 count], c2 = [p2 count]; 
    int mx = (c1>c2 ? c1 : c2), pos; 
    GMPInt *zero = [GMPInt zeroObj]; 
 
    NSMutableArray *sum = 
        [NSMutableArray arrayWithCapacity:mx]; 
 
    for(pos=0; pos<mx; pos++){

Polynomials of different degrees are handled by checking whether the current power is present in the polynomial and using a zero value otherwise.

 
        GMPInt 
            *n1 = (pos < c1 ? [p1 objectAtIndex:pos] : zero), 
            *n2 = (pos < c2 ? [p2 objectAtIndex:pos] : zero), 
            *res = [n1 add:n2]; 
        [sum addObject:res]; 
    } 
 
    return sum; 
}

The product of two polynomials is computed by iterating over all possible pairs of terms. The terms azn and bzm yield abzm+n. We start by computing the degree of the product. The coefficient array must have one more slot than this value (for the constant term). We initialize the product to be empty.

 
NSMutableArray *PolyProd(NSMutableArray *f1, NSMutableArray *f2) 
{ 
    int c1 = [f1 count], c2 = [f2 count]; 
    int cap = c1+c2-1, p, p1, p2; 
    GMPInt *zero = [GMPInt zeroObj]; 
 
    NSMutableArray *prod = 
        [NSMutableArray arrayWithCapacity:cap]; 
 
    for(p=0; p<cap; p++){ 
        [prod addObject:zero]; 
    }

We iterate over the terms of the two polynomials. We have a term for the product when both terms have non-zero coefficients.

 
    for(p1=0; p1<c1; p1++){ 
        for(p2=0; p2<c2; p2++){ 
            GMPInt 
                *n1 = [f1 objectAtIndex:p1], 
                *n2 = [f2 objectAtIndex:p2]; 
            if([n1 isZero]==NO && [n2 isZero]==NO){

If this is the case, then we compute the degree of the product and its coefficient, read the old value for the degree, increment it, and write it back into the coefficient array of the product. We return the result when we are done.

 
                int ex = p1+p2; 
                GMPInt 
                     *prev = [prod objectAtIndex:ex], 
                     *p = [n1 mul:n2], 
                     *cur = [prev add:p]; 
                [prod replaceObjectAtIndex:ex 
                       withObject:cur]; 
            } 
        } 
    } 
 
    return prod; 
}

We need to be able to compute powers of polynomials p because a term like akm in the cycle index requires us to compute (1 + zk)m. This particular example can be computed with the binomial theorem, but we use a generic routine that uses log m multiplications. The base cases are first: one when the exponent is zero and p when the exponent is one.

 
NSMutableArray *PolyPow(NSMutableArray *p, unsigned int ex) 
{ 
    NSMutableArray *factor = p, 
        *res = 
        [NSMutableArray 
            arrayWithObject:[GMPInt oneObj]]; 
 
    if(!ex){ 
        return res; 
    } 
 
    if(ex==1){ 
        return p; 
    }

We now process the bits of the exponent in turn, starting with the least significant bit. This is a folklore algorithm, too. If the bit b is set, then we include p2b in the product. We obtain p2b by squaring p2b-1 .

 
    while(ex){ 
        if(ex%2){ 
            res = PolyProd(res, factor); 
        } 
        ex >>= 1; 
        factor = PolyProd(factor, factor); 
    } 
 
    return res; 
}

The next two routines in this section are very similar. We use them either to multiply the coefficients of a polynomial by some integer value, or to divide them by an integer value.

Multiplication is first. We iterate over the terms with a loop.

 
NSMutableArray *PolyScaleMul(NSMutableArray *p, GMPInt *f) 
{ 
    int pos, mx = [p count]; 
    NSMutableArray *res = 
        [NSMutableArray arrayWithCapacity:mx]; 
    GMPInt *zero = [GMPInt zeroObj]; 
 
    for(pos=0; pos<mx; pos++){

If we find a non-zero coefficient, then we multiply it by the factor f and record the new value.

 
        GMPInt 
            *prev = [p objectAtIndex:pos], 
            *cur = zero; 
        if([prev isZero]==NO){ 
            cur = [prev mul:f]; 
        } 
        [res addObject:cur]; 
    } 
 
    return res; 
}

Division works the same way. We will use these two routines to scale terms of the cycle index by their coefficients and average the substituted cycle index at the very end.

 
NSMutableArray *PolyScaleDiv(NSMutableArray *p, GMPInt *f) 
{ 
    int pos, mx = [p count]; 
    NSMutableArray *res = 
        [NSMutableArray arrayWithCapacity:mx]; 
    GMPInt *zero = [GMPInt zeroObj]; 
 
    for(pos=0; pos<mx; pos++){ 
        GMPInt 
            *prev = [p objectAtIndex:pos], 
            *cur = zero; 
        if([prev isZero]==NO){ 
            cur = [prev div:f]; 
        } 
        [res addObject:cur]; 
    } 
 
    return res; 
}

At this point we have all the support that needs to be in place for us to be able to actually substitute 1 + z into a term from the cycle index. We must replace ak by (1 + zk) in the product k=1pakjk. This is what the next routine does. It iterates over the jk and looks for nonzero values.

 
NSMutableArray *MSetToPoly(NSMutableArray *mset) 
{ 
    NSMutableArray 
        *res = 
        [NSMutableArray 
            arrayWithObject:[GMPInt oneObj]]; 
 
    int pos, mx = [mset count]; 
 
    for(pos=0; pos<mx; pos++){ 
        int e = [[mset objectAtIndex:pos] intValue]; 
        if(e){

If it finds such a value, then it includes the term (1 + zk)jk in the product. The last step is to return the result. E.g. this will transform a1 a23 into z7 + z6 + 3z5 + 3z4 + 3z3 + 3z2 + z + 1.

 
            res = PolyProd(res, PolyPow(PolyK(pos+1), e)); 
        } 
    } 
 
    return res; 
}

We promised extensive debug facilities and hence we need to be able to convert polynomials into strings that we can output during debugging. We do this in the next routine by computing the individual terms and joining them by interpolating plus signs. We use the LATEX format, which is readable and ready for inclusion into LATEX documents.

We iterate over the terms and skip empty ones, processing the others and storing the coefficient string in the variable cfStr.

 
NSString *PolyToString(NSMutableArray *p) 
{ 
    int pos, mx = [p count]; 
    NSMutableArray *terms = 
        [NSMutableArray arrayWithCapacity:mx]; 
 
    for(pos=0; pos<mx; pos++){ 
        GMPInt *coeff = [p objectAtIndex:pos]; 
 
        if([coeff isZero]==YES){ 
            continue; 
        } 
 
        NSString 
            *cfstr = [coeff stringValue], *term;

The constant term is identical to its coefficient.

 
        if(!pos){ 
            term = cfstr; 
        }

The term for z is either z when the coefficient q is one, and qz otherwise.

 
        else if(pos==1){ 
            if([coeff isOne]==YES){ 
                term = @”z”; 
            } 
            else{ 
                term = 
                     [NSString stringWithFormat: 
                                   @”%@_z”, cfstr]; 
            } 
        }

The term for qzk where k > 1 is zk when q is one and qzk otherwise.

 
        else{ 
            if([coeff isOne]==YES){ 
                term = 
                     [NSString stringWithFormat: 
                                   @”zˆ{%d}, pos]; 
            } 
            else{ 
                term = 
                     [NSString stringWithFormat: 
                                   @”%@_zˆ{%d}, cfstr, pos]; 
            } 
        } 
 
        [terms addObject:term]; 
    }

The routine returns the terms, joined by a plus sign.

 
    return [terms componentsJoinedByString:@”_+_”]; 
}

1.5.9 Implementation III: Pólya’s theorem

Counting graphs is easy with the set of routines that we have presented and we can concentrate on the high-level aspects of the algorithm. Our program needs a bit of an interface; it will process the command line and accept two arguments, the first of which is an optional switch -d, which turns on debugging, and the second or third gives n, the number of vertices of the graphs.

We need an autorelease pool as explained earlier. We also need the command line arguments if we want to process them.

 
int main(int argc, char** argv, char **env) 
{ 
    NSAutoreleasePool *pool = [NSAutoreleasePool new]; 
 
    NSProcessInfo *procInfo = [NSProcessInfo processInfo]; 
    NSArray *args = [procInfo arguments]; 
    int ac = [args count]; 
 
    NSString *nstr;

The first thing we do before any other processing is to clear the debug flag and declare an array of classes whose memory allocation we wish to debug. In the case of mutable arrays and strings we obtain the class from an instance rather than the class itself. This is because these two classes produce instances whose class is not the same as the class that was asked to instantiate itself. We also initialize an array of allocation counters that allow us to test when allocation has doubled and empty the autorelease pool if this is the case. The variable cl is used to iterate over the array of classes.

 
    BOOL debug = NO; 
 
    Class todebug[] = { 
        [GMPInt class], 
        [[NSMutableArray array] class], 
        [[NSString stringWithFormat:@”instance_%d”, 0] class], 
        nil 
        }, *current; 
    int allocs[] = { 
        0, 0, 0 
    }; 
    GSDebugAllocationActive(YES); 
    BOOL gc; id cl; int cind;

We still have some basic command line processing to do. We require two or three arguments; everything else is an error. (The arguments include the program name.)

 
    if(ac!=2 && ac!=3){ 
        [NSException raise:NSInvalidArgumentException 
                      format:@”usage:_%@_[-d]_<n>, 
                      [procInfo processName]]; 
    }

If we have three arguments, then the second one must be the debug flag. We raise an exception if it isn’t.

 
    if(ac==3 && 
       [[args objectAtIndex:1] isEqualToString:@-d”]==NO){ 
        [NSException raise:NSInvalidArgumentException 
                      format:@”’-d’_expected,_got_%@”, 
                      [args objectAtIndex:1]]; 
    }

We turn debugging on if we did get a debug flag and extract the string for n from the third position.

 
    else if(ac==3){ 
        nstr = [args objectAtIndex:2]; 
        debug = YES; 
    }

The string for n is in the second position if there was no debug flag.

 
    else{ 
        nstr = [args objectAtIndex:1]; 
    }

We instantiate a scanner so that we can extract the value from the argument string for n. We raise an exception if the scanner could not read an integer or if it wasn’t positive. Otherwise we initialize t, which is the number of edges.

 
    NSScanner *scn = [NSScanner scannerWithString:nstr]; 
    int n = -1; 
    if([scn scanInt:&n]==NO ∣∣ n<1){ 
        [NSException raise:NSInvalidArgumentException 
                      format:@”positive_vertex_count,_please;_” 
                      @”got_%d”, n]; 
    } 
    int t = n*(n-1)/2, s;

We can now prepare for the computation. We compute n! and store the value in f, retaining f. We need f as the numerator of the coefficient of a product term in the cycle index and to check that the coefficients of those terms really do add up to n! This is the purpose of the checksum.

 
    GMPInt *f = [GMPInt factorialUI:n], 
        *vertexCheckSum = [GMPInt zeroObj]; 
    [f retain]; [vertexCheckSum retain];

We will be iterating over partitions, turning them into vertex permutations first and then edge permutations, finally substituting them with 1 + z. The variable psum holds the sum of the terms that we have processed so far. It is a polynomial and starts out being zero.

 
    NSMutableArray *psum = 
        [NSMutableArray 
            arrayWithObject: 
                [GMPInt zeroObj]]; 
    [psum retain];

We are ready for the first step, i.e. computing the partitions. We store and retain the result so that it will not be released. We must also retain the enumerator that we use to iterate over the partitions. We may now start the iteration itself.

 
    NSMutableArray *part, *parts = Partitions(n); 
    [parts retain]; 
 
    NSEnumerator *en = [parts objectEnumerator]; 
    [en retain]; 
 
    while((part = [en nextObject])!=nil){

First we convert the partition to a vertex permutation shape, then we convert that shape in turn into an edge permutation shape. We compute the denominator of the coefficient. The numerator is n!. This yields the coefficient.

 
        NSMutableArray 
            *vmset = PartitionToMSet(part), 
            *emset = VertexMSetToEdgeMSet(vmset); 
 
        GMPInt *denom = MSetToCoeffDenom(vmset), 
            *coeff = [f div:denom];

We add the coefficient to the checksum and do the substitution of 1 + zk for ak in the edge permutation multiset.

 
        [vertexCheckSum autorelease]; 
        vertexCheckSum = [vertexCheckSum add:coeff]; 
        [vertexCheckSum retain]; 
 
        NSMutableArray *poly = MSetToPoly(emset);

We temporarily leave the computation in order to output debugging information if the user did activate the debug feature. We convert the current cycle multisets for the vertex and the edge permutation to strings and output them together with the coefficient and the polynomial we obtained.

 
        if(debug==YES){ 
            NSString 
                *vprod = MSetToProduct(vmset), 
                *eprod = MSetToProduct(emset); 
 
            fprintf(stderr, ”%s_[_%s_]_[_%s_]_%s\n”, 
                     [[coeff description] cString], 
                     [vprod cString], 
                     [eprod cString], 
                     [PolyToString(poly) cString]); 
        }

Now return to the computation. We release the old sum of the substituted cycle index terms, compute the new sum, and retain it.

 
        [psum autorelease]; 
        psum = PolySum(psum, PolyScaleMul(poly, coeff)); 
        [psum retain];

Two important consistency checks are next. The cycle lengths for the vertex permutation must add up to n, anything else is an error.

 
        if((s=SumIt(vmset))!=n){ 
            [NSException 
                raise:NSInternalInconsistencyException 
                format:@”vertices_don’t_add_up:_%d,_%d”, 
                n, s]; 
        }

Similarly, the cycle lengths for the edge permutation must add up to t = 12n(n - 1).

 
        if((s=SumIt(emset))!=t){ 
            [NSException 
                raise:NSInternalInconsistencyException 
                format: @”edges_don’t_add_up:_%d,_%d”, 
                t, s]; 
        }

We leave the computation for the second time, this time in order to release objects that are no longer used. We start by setting the release flag to off and compare the number of allocated objects of the three classes that we track to the number of allocations that we have recorded for them. If this value has more than doubled, then it is time to release the pool.

 
        gc = NO; 
        current=todebug; 
        while(*current!=nil){ 
            id cl=*current++; 
            cind = current-todebug; 
 
            if(2*allocs[cind]<GSDebugAllocationCount(cl)){ 
                gc = YES; 
                break; 
            } 
        }

We release the pool if necessary and record the new values for the number of allocated objects from each class. This ends the body of the loop and hence the actual computation.

 
        if(gc==YES){ 
            if(debug==YES){ 
                fputs(”GC_...\n”, stderr); 
            } 
            [pool release]; 
            pool = [NSAutoreleasePool new]; 
            while(*current!=nil){ 
                id cl=*current++; 
                cind = current-todebug; 
 
                allocs[cind] = GSDebugAllocationCount(cl); 
            } 
        } 
    }

We subtract the actual total for the coefficients from n! and record all three values (n!, checksum, difference) in a string.

 
    GMPInt *delta = [f sub:vertexCheckSum]; 
    NSString *csinfo = 
        [NSString stringWithFormat:@”expected_%@,_” 
                   @”got_%@,_difference_%@”, 
                   [f stringValue], 
                   vertexCheckSum, delta];

We raise an exception if the checksum and n! do not match.

 
    if([delta isZero]==NO){ 
            [NSException 
                raise:NSInternalInconsistencyException 
                format: @”checksum:_%@”, csinfo]; 
    }

We print some statistics if the user requested them, namely the checksum summary and the allocation counts of all classes that we are tracking.

 
    if(debug==YES){ 
        fprintf(stderr, \nchecksum:_%s\n\n”, [csinfo cString]); 
 
        current=todebug; 
        while(*current!=nil){ 
            cl=*current++; 
            fprintf(stderr, ”%s_count_%u_total_%u_peak_%u\n”, 
                     [NSStringFromClass(cl) cString], 
                     GSDebugAllocationCount(cl), 
                     GSDebugAllocationTotal(cl), 
                     GSDebugAllocationPeak(cl)); 
        } 
    }

We are done. It remains to average the cycle index by dividing the sum of the polynomial terms by n! and print the result. This is the generating function that we are looking for. The term qzk where 1 k 1
2n(n- 1) tells us that there are q nonisomorphic graphs with k edges on n vertices.

 
    printf(”%s\n”, 
           [PolyToString(PolyScaleDiv(psum, f)) cString]); 
 
    [pool release]; 
    exit(0); 
}

1.5.10 Results

These are the first few generating functions gn for small n.

g2  =  1+ z
g3  =  1+ z + z2 + z3
g   =  1+ z + 2z2 + 3z3 + 2z4 + z5 + z6
 4              2    3    4    5    6
g5  =  1+ z + 2z + 4z + 6z + 6z + 6z
    +  4z7 + 2z8 + z9 + z10
g6  =  1+ z + 2z2 + 5z3 + 9z4 + 15z5 + 21z6
    +  24z7 + 24z8 + 21z9 + 15z10 + 9z11
         12    13   14    15
    +  5z  + 2z  + z  + z
g7  =  1+ z + 2z2 + 5z3 + 10z4 + 21z5
    +  41z6 + 65z7 + 97z8 + 131z9 + 148z10
    +  148z11 + 131z12 + 97z13 + 65z14 + 41z15
          16     17     18    19   20   21
    +  21z  + 10z  + 5z  + 2z  + z  + z
g8  =  1+ z + 2z2 + 5z3 + 11z4 + 24z5
    +  56z6 + 115z7 + 221z8 + 402z9 + 663z10
    +  980z11 + 1312z12 + 1557z13 + 1646z14
            15       16      17      18
    +  1557z  + 1312z  +980z  + 663z
    +  402z19 + 221z20 + 115z21 + 56z22
    +  24z23 + 11z24 + 5z25 + 2z26 + z27 + z28

And here is an excerpt from the result for n = 32.

...
+   34761657216148743448344973243138057890667300337466944z72
                                                        73
+   200461626459336565845980681588043820267288760177587840z
+   1138648479398347554889254519273951546678680253899786257z74
+   6371020632922419133637914798015058708350896416489205295z75
+   35117687386619298252758128999284099659276881466816638923z76
                                                           77
+   190712745660577653411399573067088056308497057681671998323z
+   1020497580980184635651931255165238794997244766608114575302z78
+   5381036625283742958677377853029609055260579018816386455844z79
+   27963157472855646088638121146117229104449204522412651394430z80

...