Pólya Enumeration

Why Pólya Enumeration Matters

A circular pattern has four positions, each colored red or blue. Rotations identify patterns, while reflections and exchanging the color names are not additional equivalences. There are six rotationally distinct patterns in total, but a designer may ask a more specific question: how many use exactly two red positions? The answer is two: red positions can be adjacent or opposite.

Burnside's lemma counts symmetry classes by averaging fixed configurations. Pólya enumeration retains the cycle lengths of each symmetry and uses generating functions to record color inventories. It can return a whole distribution instead of one total. The same idea supports catalogs of symmetric tiles, circular sequences, and substitution patterns on a fixed molecular skeleton, provided the allowed symmetries are modeled correctly.

This lesson develops the general theorem and then implements a focused case: the number of binary necklaces with each possible red-bead count, under rotations only. You will trace a cycle-index substitution, build its coefficients with dynamic programming, and compare independently runnable Python, C++, and Java programs. The polynomial recurrence connects directly to basic DP; the group action specifies which configurations count as the same object.

Core Concepts and Terminology

Let $G$ be a finite group permuting $n$ positions. It acts on colorings by moving positions while preserving color names. An orbit is one equivalence class of colorings under those transformations. For a permutation $g\in G$, let $m_\ell(g)$ be its number of cycles of length $\ell$, including cycles of length one. Every position belongs to exactly one cycle, so

$$ \sum_{\ell=1}^{n}\ell\,m_\ell(g)=n. $$

Introduce formal variables $a_1,\ldots,a_n$. The cycle index is

$$ Z_G(a_1,\ldots,a_n)=\frac{1}{\lvert G\rvert} \sum_{g\in G}\prod_{\ell=1}^{n}a_\ell^{m_\ell(g)}. $$

The subscript records cycle length; the exponent records how many such cycles occur. For rotations of four positions, the identity has four one-cycles, the half-turn has two two-cycles, and each quarter-turn has one four-cycle. Therefore

$$ Z_{C_4}=\frac{a_1^4+a_2^2+2a_4}{4}. $$

A coloring fixed by a permutation is constant on every cycle. If the available color weights are $w_1,\ldots,w_c$, a monochromatic cycle of length $\ell$ contributes one of $w_1^\ell,\ldots,w_c^\ell$. Pólya's inventory theorem makes the substitution

$$ a_\ell\longleftarrow\sum_{i=1}^{c}w_i^\ell. $$

This is the cycle-index and inventory framework described in Stanford's cycle-index notes and inventory theorem notes. Setting every color weight to one recovers the total orbit count; keeping variables distinguishes inventories.

For binary colorings, give red weight $z$ and blue weight $1$. A cycle contributes $1+z^\ell$, and the inventory polynomial is

$$ I_G(z)=Z_G(1+z,1+z^2,\ldots,1+z^n). $$

The coefficient $[z^q]I_G(z)$ counts orbits with exactly $q$ red positions. Although the cycle index has rational coefficients, the resulting inventory has nonnegative integer coefficients. Every coloring in an orbit has the same inventory, so the weight is well-defined for that orbit.

How Pólya Enumeration Works, Step by Step

The general workflow is to choose the group, classify permutation cycles, substitute color weights, combine polynomials with multiplicities, and average over the group. We specialize this workflow to rotations of $n\ge1$ positions and two named colors, without any adjacency constraint.

A rotation by $r$ positions has $d=\gcd(n,r)$ cycles, all of length $L=n/d$. Its fixed-coloring polynomial is

$$ (1+z^L)^d=\sum_{j=0}^{d}\binom{d}{j}z^{jL}. $$

Selecting $j$ cycles to be red colors $jL$ positions red. Group rotations by $L$ and let $h_L$ be their multiplicity. Then

$$ I_{C_n}(z)=\frac1n\sum_L h_L(1+z^L)^{n/L}. $$

The code works modulo $p=1{,}000{,}000{,}007$, with $1\le n\lt p$. This makes $n$ invertible modulo $p$. It computes the numerator's coefficient array first and multiplies every entry by $n^{p-2}\bmod p$ afterward.

1. Visit shifts $r=0,\ldots,n-1$. Compute $L=n/\gcd(n,r)$ and increment types[L]. The multiplicities must sum to $n$.
2. Initialize the numerator array total with $n+1$ zeros.
3. Process the distinct lengths $L$ in ascending order. Put $d=n/L$ and build the coefficients of $(1+t)^d$ in a compressed array choose of length $d+1$.
4. Begin with choose[0] = 1. For each new cycle, update indices from the newly reachable highest index down to one using choose[j] += choose[j - 1], modulo $p$.
5. For every $j=0,\ldots,d$, add types[L] * choose[j] modulo $p$ into total[j * L]. Only multiples of $L$ can occur for this cycle type.
6. Multiply all entries of total by the modular inverse of $n$ and return the array.

The compressed variable $t$ tracks red cycles; expansion replaces $t$ by $z^L$, which tracks red positions. The descending sweep is essential here: each cycle is assigned a color exactly once. This is a zero-or-one choice per cycle, unlike a DP in which an item may be reused without limit.

flowchart TD
    A["Read necklace length n"] --> B{"1 at most n and n less than p?"}
    B -- No --> C["Reject input"]
    B -- Yes --> D["Group every rotation by cycle length"]
    D --> E["Initialize numerator array to zero"]
    E --> F{"Another cycle length L?"}
    F -- Yes --> G["Set d = n divided by L; initialize choose"]
    G --> H{"Another cycle to include?"}
    H -- Yes --> I["Update choose in descending order modulo p"]
    I --> H
    H -- No --> J["Add weighted coefficients at indices j times L"]
    J --> F
    F -- No --> K["Compute inverse of n modulo p"]
    K --> M["Scale every coefficient and return inventory"]

Worked Example with a Full Trace

Use four positions, with red weight $z$ and blue weight $1$. Scan the rotations in order. Shift $0$ gives $L=1$, so the histogram becomes {1: 1}. Shift $1$ gives $L=4$, producing {1: 1, 4: 1}. Shift $2$ gives $L=2$, producing {1: 1, 4: 1, 2: 1}. Shift $3$ gives $L=4$ again, producing {1: 1, 4: 2, 2: 1}. The multiplicities sum to $4$.

Start the numerator at [0, 0, 0, 0, 0]. Each index is a red-position count.

For length $L=1$, there are $d=4$ cycles and multiplicity $1$. The compressed array begins [1, 0, 0, 0, 0]. Its four completed passes are [1, 1, 0, 0, 0], [1, 2, 1, 0, 0], [1, 3, 3, 1, 0], and [1, 4, 6, 4, 1]. For a complete inner-loop trace, the writes (index, new value) in these passes are $(1,1)$; then $(2,1),(1,2)$; then $(3,1),(2,3),(1,3)$; then $(4,1),(3,4),(2,6),(1,4)$. Every other entry retains its preceding value until its own write. Expansion changes no indices because $L=1$. The running numerator becomes [1, 4, 6, 4, 1].

For length $L=2$, there are $d=2$ cycles and multiplicity $1$. The compressed states are [1, 0, 0], [1, 1, 0], and [1, 2, 1]. The individual writes are $(1,1)$ in the first pass, then $(2,1),(1,2)$ in the second. Place these entries at red counts $0,2,4$, giving [1, 0, 2, 0, 1]. The running numerator becomes [2, 4, 8, 4, 2].

For length $L=4$, there is $d=1$ cycle and multiplicity $2$. The compressed states are [1, 0] and [1, 1]. Expansion gives [1, 0, 0, 0, 1], and the multiplicity doubles it to [2, 0, 0, 0, 2]. The final numerator is [4, 4, 8, 4, 4].

Dividing each coefficient by $4$ gives [1, 1, 2, 1, 1], or

$$ I_{C_4}(z)=\frac{(1+z)^4+(1+z^2)^2+2(1+z^4)}4 =1+z+2z^2+z^3+z^4. $$

The coefficient of $z^2$ is $2$: representatives are 0011 and 0101, with 1 denoting red. The coefficients sum to $I_{C_4}(1)=6$, agreeing with the total binary necklace count. Under the modulus, the program uses $4^{-1}\bmod p=250{,}000{,}002$; all these small results are also their exact counts.

Implementation in Python, C++, and Java

Each program returns the complete binary necklace inventory, modulo the same explicit prime. Entry $q$ counts rotational equivalence classes containing exactly $q$ red positions. Blue and red remain different named colors; exchanging them maps one inventory to its complementary inventory rather than identifying the two.

The input contract is $1\le n\lt p$, and the required arrays and computation must fit available resources. Length zero is rejected. This implementation constructs the rotation group's cycle types itself; it does not accept arbitrary groups or include reflections. The general cycle-index substitution above explains how to extend it.

All three programs print modulus = 1000000007 followed by I(1) = [1, 1], I(2) = [1, 1, 1], I(4) = [1, 1, 2, 1, 1], I(6) = [1, 1, 3, 4, 3, 1, 1], and I(8) = [1, 1, 4, 7, 10, 7, 4, 1, 1], one result per line. The notation I(n) in this output labels the coefficient array for length $n$.

Python uses exact integers, C++ uses long long, and Java uses BigInteger. C++ reduces each weighted product before adding it. Each product is at most $(p-1)^2=1{,}000{,}000{,}012{,}000{,}000{,}036$, safely below the signed $64$-bit limit, and each subsequent sum is at most $2p-2$. The programs never use floating-point division. In Java, save the complete panel as Main.java.

from math import gcd

MOD = 1_000_000_007


def necklace_inventory(n: int) -> list[int]:
    if not 1 <= n < MOD:
        raise ValueError("require 1 <= n < MOD")
    types = {}
    for shift in range(n):
        length = n // gcd(n, shift)
        types[length] = types.get(length, 0) + 1
    total = [0] * (n + 1)
    for length in sorted(types):
        cycles = n // length
        choose = [0] * (cycles + 1)
        choose[0] = 1
        for processed in range(cycles):
            for red in range(processed + 1, 0, -1):
                choose[red] = (choose[red] + choose[red - 1]) % MOD
        for red_cycles in range(cycles + 1):
            index = red_cycles * length
            weighted = types[length] * choose[red_cycles] % MOD
            total[index] = (total[index] + weighted) % MOD
    inverse_n = pow(n, MOD - 2, MOD)
    return [value * inverse_n % MOD for value in total]


def main() -> None:
    print(f"modulus = {MOD}")
    for n in [1, 2, 4, 6, 8]:
        print(f"I({n}) = {necklace_inventory(n)}")


if __name__ == "__main__":
    main()

#include <iostream>
#include <map>
#include <stdexcept>
#include <vector>

constexpr long long MOD = 1000000007LL;

int gcd_nonnegative(int a, int b) {
    while (b != 0) {
        int remainder = a % b;
        a = b;
        b = remainder;
    }
    return a;
}

long long mod_pow(long long base, long long exponent) {
    long long result = 1;
    while (exponent > 0) {
        if (exponent % 2 == 1) result = result * base % MOD;
        base = base * base % MOD;
        exponent /= 2;
    }
    return result;
}

std::vector<long long> necklace_inventory(int n) {
    if (n < 1 || n >= MOD) {
        throw std::invalid_argument("require 1 <= n < MOD");
    }
    std::map<int, int> types;
    for (int shift = 0; shift < n; ++shift) {
        int length = n / gcd_nonnegative(n, shift);
        ++types[length];
    }
    std::vector<long long> total(n + 1, 0);
    for (const auto& type : types) {
        int length = type.first;
        int cycles = n / length;
        std::vector<long long> choose(cycles + 1, 0);
        choose[0] = 1;
        for (int processed = 0; processed < cycles; ++processed) {
            for (int red = processed + 1; red >= 1; --red) {
                choose[red] = (choose[red] + choose[red - 1]) % MOD;
            }
        }
        for (int red_cycles = 0; red_cycles <= cycles; ++red_cycles) {
            int index = red_cycles * length;
            long long weighted = type.second * choose[red_cycles] % MOD;
            total[index] = (total[index] + weighted) % MOD;
        }
    }
    long long inverse_n = mod_pow(n, MOD - 2);
    for (long long& value : total) value = value * inverse_n % MOD;
    return total;
}

int main() {
    std::cout << "modulus = " << MOD << '\n';
    for (int n : {1, 2, 4, 6, 8}) {
        const auto inventory = necklace_inventory(n);
        std::cout << "I(" << n << ") = [";
        for (std::size_t i = 0; i < inventory.size(); ++i) {
            if (i != 0) std::cout << ", ";
            std::cout << inventory[i];
        }
        std::cout << "]\n";
    }
    return 0;
}

import java.math.BigInteger;
import java.util.Arrays;
import java.util.Map;
import java.util.TreeMap;

public class Main {
    private static final int PRIME = 1_000_000_007;
    private static final BigInteger MOD = BigInteger.valueOf(PRIME);

    private static int gcdNonnegative(int a, int b) {
        while (b != 0) {
            int remainder = a % b;
            a = b;
            b = remainder;
        }
        return a;
    }

    public static BigInteger[] necklaceInventory(int n) {
        if (n < 1 || n >= PRIME) {
            throw new IllegalArgumentException("require 1 <= n < MOD");
        }
        Map<Integer, Integer> types = new TreeMap<>();
        for (int shift = 0; shift < n; ++shift) {
            int length = n / gcdNonnegative(n, shift);
            types.merge(length, 1, Integer::sum);
        }
        BigInteger[] total = new BigInteger[n + 1];
        Arrays.fill(total, BigInteger.ZERO);
        for (Map.Entry<Integer, Integer> type : types.entrySet()) {
            int length = type.getKey();
            int cycles = n / length;
            BigInteger[] choose = new BigInteger[cycles + 1];
            Arrays.fill(choose, BigInteger.ZERO);
            choose[0] = BigInteger.ONE;
            for (int processed = 0; processed < cycles; ++processed) {
                for (int red = processed + 1; red >= 1; --red) {
                    choose[red] = choose[red].add(choose[red - 1]).mod(MOD);
                }
            }
            BigInteger multiplicity = BigInteger.valueOf(type.getValue());
            for (int redCycles = 0; redCycles <= cycles; ++redCycles) {
                int index = redCycles * length;
                BigInteger weighted = multiplicity.multiply(choose[redCycles]).mod(MOD);
                total[index] = total[index].add(weighted).mod(MOD);
            }
        }
        BigInteger inverseN = BigInteger.valueOf(n).modPow(
                MOD.subtract(BigInteger.valueOf(2)), MOD);
        for (int i = 0; i <= n; ++i) {
            total[i] = total[i].multiply(inverseN).mod(MOD);
        }
        return total;
    }

    public static void main(String[] args) {
        System.out.println("modulus = " + MOD);
        for (int n : new int[]{1, 2, 4, 6, 8}) {
            System.out.println("I(" + n + ") = " + Arrays.toString(necklaceInventory(n)));
        }
    }
}

Complexity, Correctness, and Invariants

For this implementation, let $\tau(n)$ be the number of divisors of $n$. Rotation classification uses $n$ greatest-common-divisor computations and dictionary operations, taking $O(n\log(n+1))$ time under ordinary integer and map bounds. Every occurring cycle length divides $n$.

A type with $d$ cycles performs exactly $d(d+1)/2$ coefficient updates. Across all types, the dominant work is $\Theta(\sum_{d\mid n}d^2)$. This sum is at least $n^2$, because the identity contributes $d=n$, and at most

$$ n^2\sum_{L=1}^{\infty}\frac1{L^2}\lt2n^2. $$

Thus the DP work is $\Theta(n^2)$. The total bound is $O(n^2+\log p)$ time, including the inverse, and $\Theta(n)$ space for the output, temporary coefficients, and type map. With the fixed prime, the time is $\Theta(n^2)$. There is no configuration-dependent branching: for a fixed length, best, average, and worst cases have the same operation count. An average across different lengths would require a distribution. These are arithmetic-operation bounds for bounded modular values, not arbitrary-size integer bit bounds.

Correctness has two levels. A coloring fixed by a permutation assigns one color to each complete cycle, giving the substituted product. Averaging these fixed-coloring polynomials counts each orbit once because its weight is constant across the orbit: the weighted version of Burnside's double-counting proof applies coefficient by coefficient. For the implementation, the compressed recurrence produces the binomial coefficients, expansion maps red cycles to red positions, and multiplicities restore all group elements. Multiplication by the valid inverse performs the final averaging modulo $p$.

Variants, Pitfalls, and Optimizations

More colors and fixed inventories. Give three colors weights $x,y,z$. Replace $a_L$ by $x^L+y^L+z^L$ and extract $[x^r y^g z^b]$ for prescribed counts. A cycle of length $L$ contributes $L$ positions of one color, not one position. The polynomial dimension and coefficient storage can grow quickly; track only the counts needed, with one color often determined by the total.

Recover a scalar count. Set all inventory variables to one. With $c$ unrestricted colors, each cycle-index variable becomes $c$, yielding $Z_G(c,\ldots,c)$. For the binary program, summing the returned coefficients modulo $p$ recovers the binary necklace total. Retaining only that total would discard the inventory information this lesson computes.

Choose the correct group. Include reflections to count bracelets; include other position permutations only when the application declares them equivalent. Four rotations form $C_4$, while a square's rotations and reflections form $D_4$ with eight elements. Renaming colors is another equivalence and is not automatically included by the ordinary position-action substitution. For constrained colorings such as proper graph colorings, plain independent cycle choices may fail; the constraints must be preserved and counted within each fixed set.

Build cyclic types directly. The cycle index of the rotation group also has the divisor form

$$ Z_{C_n}=\frac1n\sum_{L\mid n}\varphi(L)a_L^{n/L}, $$

where $\varphi$ is Euler's totient function. Thus $h_L=\varphi(L)$. Factoring $n$ and generating its divisors can replace the scan of all shifts. Include those preparation costs when comparing implementations.

Generate binomial coefficients faster. The teaching code builds each row by additions to expose the DP invariant. Since $n\lt p$, factorial and inverse-factorial tables can generate $\binom{d}{j}\bmod p$ in constant time per entry after preprocessing. A multiplicative recurrence with valid modular inverses is another option. Neither approach permits ordinary integer division of already-reduced residues.

Keep cycle lengths in the substitution. The tempting replacement $a_2\leftarrow1+z$ counts red two-cycles rather than red positions. The correct replacement is $1+z^2$. Applied to the half-turn, the incorrect product $(1+z)^2$ invents a contribution at red count one; no half-turn-invariant coloring can have only one red position.

Preserve descending updates and multiplicities. An ascending compressed sweep can reuse a cycle within one pass, producing coefficients that no longer equal the binomial row. Dropping the second quarter-turn also corrupts the numerator. Finally, if the group order is not invertible under the chosen modulus, divide exact integer coefficients before reducing, or use an appropriate larger-modulus method; this code rejects lengths outside its stated invertible range.

Practice Problems and Real-World Use

Start each problem by specifying what is permuted, which transformations are allowed, and what the variables measure. Compute the cycle index or enough of its terms to extract the target coefficient. The small targets below can be verified by canonicalizing all allowed colorings under the stated symmetries.

For an additional check, exchange red and blue. This gives a bijection between inventories with $q$ and $n-q$ red positions, so every returned array must be symmetric. Its two endpoint coefficients must both be $1$, its entries must sum to the total necklace count modulo $p$, and the exhaustive verifier must agree coefficient by coefficient for short lengths.