Pólya Enumeration Theorem

July 2026 · 28 min read

From Burnside to Pólya

In the previous post we met Burnside’s lemma: the number of distinct colourings of a set X under a group G of symmetries is

|X/G|  =  (1/|G|)  Σ_g∈G |Fix(g)|

Burnside tells us how many distinct colourings exist, but nothing about their structure. What if we want to know how many distinct necklaces use exactly 3 red beads and 2 blue beads? Burnside alone cannot answer that—it only yields a single count.

Pólya’s Enumeration Theorem (PET) extends Burnside’s lemma by introducing a weight for each colour. Instead of merely counting fixed points, we substitute colour–weight polynomials into the cycle index of the group and obtain a generating function whose coefficients tell us the number of distinct colourings by weight distribution.

The central object that makes PET work is the cycle index of the symmetry group. Once you can compute the cycle index, the rest is pure substitution. Let’s define it precisely.

Cycle Index of a Permutation Group

Definition

Let G be a permutation group acting on a set of n elements. Every permutation σ∈G can be decomposed into disjoint cycles. If σ has j_k cycles of length k (for k = 1, 2, …, n), we write its cycle type as (j₁, j₂, …, j_n). Note that Σ_k k·j_k = n.

The cycle index of G is the polynomial in indeterminates s₁, s₂, …, s_n:

Z(G; s₁, s₂, …, s_n)  =  (1/|G|)  Σ_σ∈G  s₁^j₁(σ)  s₂^j₂(σ)  …  s_n^j_n(σ)

In other words, each permutation contributes a monomial that encodes its cycle structure, and the cycle index is the average of all these monomials.

Cycle index of the cyclic group C_n

The cyclic group C_n consists of the n rotations r⁰, r¹, …, rⁿ⁻¹ of n objects arranged in a circle. The rotation rⁱ has gcd(i, n) cycles, each of length n/gcd(i, n). Grouping rotations by gcd with n and applying Euler’s totient function φ gives the elegant closed form:

Z(C_n)  =  (1/n) Σ_d|n φ(d) · s_d^n/d

Derivation. Consider rotation rⁱ. Let d = gcd(i, n). This permutation consists of exactly d cycles, each of length n/d, so it contributes s_n/d^d. For a given divisor d of n, how many values of i (0 ≤ i < n) satisfy gcd(i, n) = d? There are φ(n/d) such values. Substituting d’ = n/d and relabelling:

Z(C_n) = (1/n) Σ_d|n φ(n/d) · s_n/d^d = (1/n) Σ_d’|n φ(d’) · s_d’^n/d’

Cycle index of the dihedral group D_n

The dihedral group D_n adds n reflections to C_n, giving |D_n| = 2n. The reflection cycle types depend on parity:

So the full cycle index of D_n is:

Z(D_n)  =  ½ Z(C_n)  +  (reflection terms) / (2n)

Worked example: Z(C₄) and Z(D₄)

For n = 4 the divisors are 1, 2, 4.

Rotations (C₄):

• r⁰ (identity): 4 fixed points → s₁⁴
• r¹: one 4-cycle → s₄
• r²: two 2-cycles → s₂²
• r³: one 4-cycle → s₄

Z(C₄) = (1/4)(s₁⁴ + s₂² + 2s₄)

Using the φ-formula: φ(1)·s₁⁴ + φ(2)·s₂² + φ(4)·s₄¹ = 1·s₁⁴ + 1·s₂² + 2·s₄ ✔

Reflections (n = 4, even):

• 2 reflections through vertices: each contributes s₁² s₂
• 2 reflections through edges: each contributes s₂²

Z(D₄) = (1/8)(s₁⁴ + 2s₄ + 3s₂² + 2s₁² s₂)

The Pólya Enumeration Theorem

Unweighted version (counting distinct colourings)

Suppose we want to colour n positions using k colours under a symmetry group G. The total number of distinct colourings is obtained by substituting s_i = k for all i:

|colourings|  =  Z(G; k, k, …, k)

Weighted (inventory) version

Now let each colour c have a formal weight w(c). Define the figure-counting series:

f(x)  =  Σ_c w(c)

The Pólya Enumeration Theorem states that the pattern inventory (generating function counting distinct colourings by weight) is:

PET(G)  =  Z(G; s₁→f(x), s₂→f(x²), …, s_n→f(xⁿ))

More explicitly, for each indeterminate s_i, substitute:

s_i  →  Σ_c w(c)ⁱ

Example. Two colours {R, B} with weights w(R) = r, w(B) = b. Then s_i → rⁱ + bⁱ. The resulting polynomial in r and b enumerates distinct colourings: the coefficient of r^ab^c is the number of distinct colourings using exactly a reds and c blues.

Connection to Burnside

Setting every w(c) = 1 collapses the inventory to a single number. Each s_i becomes 1ⁱ + 1ⁱ + … = k (for k colours), and PET reduces to the unweighted formula, which is precisely Burnside’s lemma. So PET is a strict generalisation.

Proof sketch

The proof follows the same orbit-counting idea as Burnside. A colouring f: X → C is fixed by σ iff f is constant on every cycle of σ. The weight of such a colouring is the product of w(c_cycle)^|cycle| over all cycles. Summing these weights over all fixed colourings and averaging over G yields the substitution formula above.

Necklace and Bracelet Counting

The most classic application of Pólya’s theorem in competitive programming is counting necklaces (rotational equivalence) and bracelets (rotational + reflectional equivalence).

Necklaces — cyclic group C_n

A necklace is an equivalence class of colourings under rotation. Using the cycle index of C_n with s_d = k:

N(n, k)  =  (1/n) Σ_d|n φ(d) · k^n/d

Example: n = 6, k = 2.

Divisors of 6: 1, 2, 3, 6.

• φ(1)·2⁶ = 1 · 64 = 64
• φ(2)·2³ = 1 · 8 = 8
• φ(3)·2² = 2 · 4 = 8
• φ(6)·2¹ = 2 · 2 = 4

N(6, 2) = (64 + 8 + 8 + 4) / 6 = 84 / 6 = 14 distinct necklaces.

Bracelets — dihedral group D_n

A bracelet is an equivalence class under both rotation and reflection (you can flip the necklace over). We add the reflection contributions:

B(n, k)  =  ½ N(n, k)  +  { (k^⌈n/2⌉)/2 if n is odd;  (k^n/2+1 + k^n/2)/4 if n is even }

Equivalently, for n odd:

B(n, k) = ½ N(n, k) + ½ k^(n+1)/2

For n even:

B(n, k) = ½ N(n, k) + ¼ (k^{n/2 + 1} + k^n/2)

Example: n = 6, k = 2.

N(6, 2) = 14 (from above). n is even, so reflection term = (2⁴ + 2³)/4 = (16 + 8)/4 = 6.

B(6, 2) = 14/2 + 6 = 7 + 6 = 13.

Wait—let’s double-check with the explicit D₆ cycle index.

Verification for n = 6, k = 2:

D₆ has 12 elements: 6 rotations + 6 reflections.

Rotation contributions (sum of k^cycles): 2⁶ + 2¹ + 2² + 2³ + 2² + 2¹ = 64 + 2 + 4 + 8 + 4 + 2 = 84.

Reflection contributions (n = 6, even): 3 reflections through vertices give 2⁴ = 16 each → 48; 3 reflections through edges give 2³ = 8 each → 24. Total reflections: 48 + 24 = 72.

B(6, 2) = (84 + 72) / 12 = 156 / 12 = 13. ✔

Necklaces with fixed colour counts

Use the weighted PET. For a necklace of n = 4 beads with 2 colours {R, B} where we want exactly 2 reds and 2 blues, substitute s_i → rⁱ + bⁱ into Z(C₄):

Z(C₄) = (1/4)(s₁⁴ + s₂² + 2s₄)

= (1/4)((r+b)⁴ + (r²+b²)² + 2(r⁴+b⁴))

The coefficient of r²b² in this expansion:

• From (r+b)⁴: C(4,2) = 6
• From (r²+b²)² = r⁴ + 2r²b² + b⁴: coefficient is 2
• From 2(r⁴+b⁴): coefficient of r²b² is 0

Total = (6 + 2 + 0)/4 = 2 necklaces with exactly 2R, 2B. They are RRBB and RBRB (up to rotation).

C++ Implementation

Let’s build a complete modular-arithmetic toolkit for necklace and bracelet counting. All computations are mod a prime p to handle the huge numbers that arise when n and k are large (up to 10⁹).

Core utilities: modular exponentiation, inverse, Euler totient

#include <bits/stdc++.h>
using namespace std;

using ll = long long;
const ll MOD = 1e9 + 7;

ll power(ll base, ll exp, ll mod = MOD) {
    ll result = 1;
    base %= mod;
    while (exp > 0) {
        if (exp & 1) result = result * base % mod;
        base = base * base % mod;
        exp >>= 1;
    }
    return result;
}

ll modinv(ll a, ll mod = MOD) {
    return power(a, mod - 2, mod); // Fermat's little theorem, mod must be prime
}

ll euler_totient(ll n) {
    ll result = n;
    for (ll p = 2; p * p <= n; p++) {
        if (n % p == 0) {
            while (n % p == 0) n /= p;
            result -= result / p;
        }
    }
    if (n > 1) result -= result / n;
    return result;
}

Counting necklaces mod p

// Returns the number of distinct necklaces with n beads
// and k colours, modulo MOD.
ll count_necklaces(ll n, ll k) {
    // N(n,k) = (1/n) * sum_{d|n} phi(d) * k^{n/d}
    ll inv_n = modinv(n);
    ll total = 0;

    // Enumerate divisors of n
    for (ll d = 1; d * d <= n; d++) {
        if (n % d != 0) continue;

        // d is a divisor
        ll phi_d = euler_totient(d);
        total = (total + phi_d % MOD * power(k, n / d)) % MOD;

        // n/d is also a divisor (if different from d)
        if (d != n / d) {
            ll d2 = n / d;
            ll phi_d2 = euler_totient(d2);
            total = (total + phi_d2 % MOD * power(k, d)) % MOD;
        }
    }

    return total % MOD * inv_n % MOD;
}

// Example usage
int main() {
    ll n, k;
    cin >> n >> k;
    cout << count_necklaces(n, k) << "\n";
    // n=6, k=2 => 14
}

Counting bracelets mod p

ll count_bracelets(ll n, ll k) {
    // B(n,k) = (1/2) * N(n,k) + reflection_term
    // where reflection_term depends on parity of n.

    ll necklaces = count_necklaces(n, k);
    ll inv2 = modinv(2);
    ll reflection;

    if (n % 2 == 1) {
        // n odd: n reflections, each with 1 fixed point
        // and (n-1)/2 transpositions => k^{(n+1)/2} per reflection
        // reflection_term = (1/2) * k^{(n+1)/2}
        reflection = inv2 % MOD * power(k, (n + 1) / 2) % MOD;
    } else {
        // n even: n/2 reflections through vertices (k^{n/2+1})
        //       + n/2 reflections through edges   (k^{n/2})
        // reflection_term = (1/4)(k^{n/2+1} + k^{n/2})
        ll inv4 = modinv(4);
        ll a = power(k, n / 2 + 1);
        ll b = power(k, n / 2);
        reflection = inv4 % MOD * ((a + b) % MOD) % MOD;
    }

    return (necklaces % MOD * inv2 % MOD + reflection) % MOD;
}

// Example: n=6, k=2 => 13

Weighted Pólya counting

When we need the full pattern inventory (e.g. how many necklaces have exactly a beads of colour 1 and b beads of colour 2), we must work with polynomials. Here is a simplified version for two colours under C_n:

// Polynomial multiplication (truncated to degree n)
using Poly = vector<ll>;

Poly poly_mult(const Poly& a, const Poly& b, int max_deg) {
    Poly c(max_deg + 1, 0);
    for (int i = 0; i < (int)a.size() && i <= max_deg; i++)
        for (int j = 0; j < (int)b.size() && i + j <= max_deg; j++)
            c[i + j] = (c[i + j] + a[i] * b[j]) % MOD;
    return c;
}

// Raise polynomial p to the e-th power, keeping terms up to max_deg
Poly poly_pow(Poly p, int e, int max_deg) {
    Poly result(max_deg + 1, 0);
    result[0] = 1;
    while (e > 0) {
        if (e & 1) result = poly_mult(result, p, max_deg);
        p = poly_mult(p, p, max_deg);
        e >>= 1;
    }
    return result;
}

// Weighted necklace count for 2 colours, n beads.
// Returns polynomial where coeff[a] = number of distinct necklaces
// with exactly a beads of colour 1 (and n-a of colour 2).
Poly weighted_necklaces(int n) {
    ll inv_n = modinv(n);
    Poly total(n + 1, 0);

    for (int d = 1; d * d <= n; d++) {
        if (n % d != 0) continue;

        auto process = [&](int div) {
            ll phi_val = euler_totient(div) % MOD;
            int cycle_len = n / div;
            int num_cycles = div;
            // s_cycle_len -> (x^cycle_len + 1) represents
            // colour1^cycle_len + colour2^cycle_len
            Poly base(cycle_len + 1, 0);
            base[0] = 1;                // colour 2 contributes x^0
            base[cycle_len] = 1;        // colour 1 contributes x^cycle_len
            Poly contrib = poly_pow(base, num_cycles, n);
            for (int i = 0; i <= n; i++)
                total[i] = (total[i] + phi_val * contrib[i]) % MOD;
        };

        process(d);
        if (d != n / d) process(n / d);
    }

    for (int i = 0; i <= n; i++)
        total[i] = total[i] * inv_n % MOD;

    return total;
}

For n = 4 the output polynomial is [1, 1, 2, 1, 1], meaning 1 necklace with 0 reds, 1 with 1 red, 2 with 2 reds, 1 with 3 reds, 1 with 4 reds. Total = 6 = N(4, 2). ✔

Generating Functions with Pólya

One of the most powerful features of PET is that it naturally produces generating functions. By assigning each colour a formal variable, the cycle index substitution yields a multivariate polynomial whose coefficients enumerate distinct structures by their “weight profile.”

Setup

Given colours C = {c₁, c₂, …, c_k} with formal weights w(c_i) = x_i, the pattern inventory is:

P(x₁, …, x_k)  =  Z(G; s_i → x₁ⁱ + x₂ⁱ + … + x_kⁱ)

The coefficient of x₁^a₁ x₂^a₂ … x_k^a_k counts the number of distinct colourings using exactly a_i copies of colour c_i.

Example: graph colouring (cube faces)

Consider colouring the 6 faces of a cube with colours {R, B, G}. The rotation group of the cube has 24 elements. Its cycle index (acting on faces) is:

Z = (1/24)(s₁⁶ + 6s₁²s₄ + 3s₁²s₂² + 8s₃² + 6s₂³)

Substituting s_i → 3 (for k = 3 colours, unweighted): Z(3, 3, 3, 3, 3, 3) = (729 + 6·9·3 + 3·9·9 + 8·9 + 6·27)/24 = (729 + 162 + 243 + 72 + 162)/24 = 1368/24 = 57.

For the weighted version (tracking colour counts), substitute s_i → rⁱ + bⁱ + gⁱ. The resulting polynomial in r, b, g gives the complete breakdown. For example, the coefficient of r³b²g tells you how many cubes have 3 red, 2 blue, and 1 green face.

Chemical isomers

PET was originally motivated by chemistry. To count distinct chemical isomers of a substituted benzene ring C₆H_6−mX_m (e.g. how many dichlorobenzenes exist), we use Z(D₆) and substitute with two “colours” (H and X). The answer for m = 2 is the coefficient of x² in the pattern inventory, which equals 3 (ortho, meta, para).

CP Problem Walkthrough

Problem: Coloured Necklace (CSES-style)

Statement. Given n (1 ≤ n ≤ 10⁹) and k (1 ≤ k ≤ 10⁹), compute the number of distinct necklaces with n beads and k colours, modulo 10⁹+7. Two necklaces are identical if one can be obtained from the other by rotation.

Analysis

This is a direct application of the necklace formula:

N(n, k) = (1/n) Σ_d|n φ(d) · k^n/d  (mod 10⁹+7)

The key algorithmic challenges are:

1. Enumerating divisors of n in O(√n) time.
2. Computing φ(d) for each divisor d. Since d can be up to 10⁹, we use the O(√d) trial-division totient.
3. Modular exponentiation k^n/d mod p in O(log(n/d)) time.
4. Modular inverse of n using Fermat’s little theorem.

The total number of divisors of n is at most about 1344 for n ≤ 10⁹ (a well-known bound), and computing φ for each costs O(√d), so the overall complexity is roughly O(√n · d(n)) which is very fast.

Full solution

#include <bits/stdc++.h>
using namespace std;
using ll = long long;
const ll MOD = 1e9 + 7;

ll power(ll b, ll e, ll m = MOD) {
    ll r = 1; b %= m;
    for (; e; e >>= 1) {
        if (e & 1) r = r * b % m;
        b = b * b % m;
    }
    return r;
}

ll modinv(ll a) { return power(a, MOD - 2); }

ll euler_phi(ll n) {
    ll res = n;
    for (ll p = 2; p * p <= n; p++) {
        if (n % p == 0) {
            while (n % p == 0) n /= p;
            res -= res / p;
        }
    }
    if (n > 1) res -= res / n;
    return res;
}

int main() {
    ios::sync_with_stdio(false);
    cin.tie(nullptr);

    ll n, k;
    cin >> n >> k;

    // Collect all divisors of n
    vector<ll> divs;
    for (ll d = 1; d * d <= n; d++) {
        if (n % d == 0) {
            divs.push_back(d);
            if (d != n / d) divs.push_back(n / d);
        }
    }

    ll ans = 0;
    for (ll d : divs) {
        ll phi = euler_phi(d) % MOD;
        ll contrib = phi * power(k, n / d) % MOD;
        ans = (ans + contrib) % MOD;
    }

    ans = ans % MOD * modinv(n) % MOD;
    cout << ans << "\n";
}

Testing

Variant: bracelet counting

If the problem instead asks for bracelets (rotation + reflection), simply add the reflection terms as shown in the count_bracelets function above. Many contest problems (e.g. “NECKLACE” on SPOJ) explicitly state whether flipping is allowed—read carefully!

Variant: colourings with forbidden patterns

Some harder problems restrict which colours may be adjacent. PET alone does not handle adjacency constraints. For those, combine PET with transfer matrix / DP approaches: use the matrix to count valid colourings fixed by each group element, then average via Burnside.

Summary

Key formulas to memorise for contests:

• Cycle index of C_n: Z(C_n) = (1/n) Σ_d|n φ(d) · s_d^n/d
• Necklaces: N(n, k) = (1/n) Σ_d|n φ(d) · k^n/d
• Bracelets (odd n): B = N/2 + k^(n+1)/2/2
• Bracelets (even n): B = N/2 + (k^n/2+1 + k^n/2)/4
• Weighted PET: substitute s_i → Σ_c w(c)ⁱ

Implementation checklist:

1. Enumerate divisors of n in O(√n)
2. Compute φ(d) for each divisor via trial division
3. Use fast modular exponentiation for k^n/d
4. Use Fermat inverse for 1/n (mod prime)
5. For bracelets, add reflection terms and divide by 2

In the next post, we’ll take a step back and cover the group theory basics that underpin Burnside and Pólya, including formal definitions of groups, orbits, and stabilisers.