∑MathAdvanced

Pólya Enumeration

Key Points

•
Pólya Enumeration Theorem generalizes Burnside’s Lemma by turning counting under symmetry into a polynomial substitution problem.
•
Every group action has a cycle index polynomial Z(G) that encodes the cycle structure of each group element.
•
To count colorings with m colors, substitute $p_{j}$ -> m in Z(G); to count weighted colorings, substitute $p_{j}$ -> $x_{1}$ ^j + $x_{2}$ ^j + ... + $x_{m}$ ^j.
•
Coefficients after substitution give counts refined by weights (e.g., the number of black beads), making PET a powerful tool for distribution queries.
•
For cyclic and dihedral groups (necklaces and bracelets), PET yields simple closed forms that are fast to compute.
•
PET integrates naturally with generating functions, enabling coefficient extraction for exact distributions.
•
Implementation hinges on computing cycle types and multiplying simple polynomials, which can be optimized using binomial expansions.
•
Common pitfalls include wrong cycle types for reflections, forgetting the 1/|G| average, and mishandling integer divisions or overflow.

Prerequisites

→Basic group theory (groups, permutations, cycles) — PET relies on understanding how group elements act as permutations and how to read cycle structures.
→Burnside's Lemma — PET generalizes Burnside’s averaging; you must know why averaging fixed colorings yields orbit counts.
→Generating functions and coefficient extraction — Weighted counting in PET is performed via polynomial substitutions and reading coefficients.
→Big integer arithmetic — Counts can be large; precise arithmetic avoids overflow and preserves exact divisibility by |G|.
→Number theory (Euler's totient, divisors) — Closed forms for cyclic/dihedral groups use divisor sums and \varphi(d).
→Polynomial algebra (convolution/binomial theorem) — Efficient PET implementations rely on multiplying sparse polynomials and using binomial expansions.

Detailed Explanation

Tap terms for definitions

01Overview

Pólya Enumeration Theorem (PET) is a method for counting distinct colorings of objects when some colorings are considered the same due to symmetry (rotations, reflections, etc.). While Burnside’s Lemma already tells us to average the number of fixed colorings over all group elements, PET packages this averaging into a single algebraic object: the cycle index polynomial Z(G). Each group element contributes a monomial reflecting how it permutes positions (its cycle structure), and Z(G) is the average of these monomials. The beauty appears when we substitute each cycle variable p_j by a color-weight expression like x_1^j + x_2^j + ... + x_m^j. After substitution, the resulting polynomial’s total value counts distinct colorings with m colors, and its coefficients answer refined questions—like how many orbits have exactly b black beads. PET is particularly effective for problems on necklaces (cyclic symmetry), bracelets (dihedral symmetry), grids or polyhedra faces under rotation groups, and any finite permutation group action. It transforms a combinatorial problem into algebraic manipulations, often reducing to simple binomial or multinomial expansions and coefficient extraction.

02Intuition & Analogies

Imagine you’re decorating seats around a circular table with colored cushions. Two arrangements that differ only by rotating the table are not really different. Burnside’s Lemma says: to count truly different decorations, for each rotation, count how many decorations look the same after that rotation, and then average. Now, instead of recounting every time, Pólya suggests a toolkit: label every rotation by how it groups seats into cycles (orbits of positions under that rotation). A rotation that sends seat 1 to seat 4, 4 to 7, and so on in a loop forms a cycle; coloring that cycle with a single color ensures the whole cycle is fixed by that rotation. If a rotation breaks the table into g cycles each of length L, then fixed colorings are obtained by choosing one color per cycle, giving m^g choices for m colors. PET packages this into variables p_1, p_2, p_3, ... where p_j stands for “a cycle of length j.” Each symmetry contributes the product p_{len(cycle)} over all its cycles. Average these and you get Z(G), the cycle index. Finally, say each color c contributes a weight x_c; then a j-cycle colored by color c contributes x_c^j. Summing over colors makes a j-cycle contribute x_1^j + x_2^j + ... + x_m^j. Substitute p_j by this sum and multiply across cycles; the result counts all orbits while keeping track of weights. If you only care about how many black cushions you used, give black weight z and every other color weight 1 (or k-1 if you aggregate all non-black colors). Then every j-cycle either adds z^j (if black) or a constant. Expanding and reading the coefficient of z^b tells you how many inequivalent arrangements use exactly b black cushions.

03Formal Definition

Let a finite group G act on a finite set X of n positions. For g

\in

G, write its action on X as a permutation and let

a_{j}

(g) be the number of cycles of length j in this permutation. The cycle index polynomial of G is defined as

\n \nZ

(G) =

\frac{1}{∣ G ∣}

\sum_{g \in G}

\prod_{j \geq 1}

j^{a_{j} (g)}

\n \nGiven

a (possibly weighted) color set C and a map w: C

\to

R (weights, often monomials in variables), the Pólya substitution replaces each

p_{j}

\n \nP_{j}

\sum_{c \in C}

w(c)^j.

\n \nThe

weighted orbit enumerator is then

\n \nF_{G}

(w) = Z(G)

\big

|_{

p_{j}

\to

P_{j}

} =

\frac{1}{∣ G ∣}

\sum_{g \in G}

\prod_{j \geq 1}

\left

(

\sum_{c \in C}

w(c)^j

\right

)^{

a_{j}

(g)}.

\n \nIn

particular, if |C| = m and unweighted (w(c) = 1), then

P_{j} = m

and the number of orbits (inequivalent colorings) is

F_{G}

\frac{1}{∣ G ∣}

\sum_{g \in G}

m^{\sum_{j} a_{j} (g)}

where

\sum_{j}

a_{j}

(g) is the number of cycles of g. For a refinement tracking the count of a distinguished color B with weight z and all other (k-1) indistinguishable colors with weight 1, use

P_{j} = z^{j}

+ (k-1). The coefficient of

z^{b}

F_{G}

is the number of orbits with exactly b positions colored B.

04When to Use

Use PET whenever you need to count colorings or labelings up to symmetry and especially when you want refined counts by color multiplicities. Classical applications include: (1) Necklaces (C_n actions) and bracelets (D_n actions): count with k colors or count distributions of a special color. (2) Colorings of faces/edges/vertices of polyhedra under their rotation groups (e.g., cube faces with 24 rotations). (3) Grid colorings under the dihedral group of the square. (4) Counting functions up to permutation of coordinates or counting multiset partitions under symmetry. (5) When you need more than just the total number—e.g., the number of inequivalent colorings with exactly b red cells—PET’s weighted substitution provides a direct, often faster, path than brute force. PET is ideal if (a) the group is small or has a closed-form cycle index (like C_n or D_n), (b) you can compute cycle types efficiently (e.g., via generators), and (c) you can afford polynomial multiplications up to degree n. It’s less appropriate if the action is enormous and no compact cycle index or generator description is available, or if constraints break the product structure (e.g., local adjacency constraints that aren’t color-symmetric).

⚠️Common Mistakes

• Forgetting the average: Burnside and PET always divide by |G|. Accumulating cycle contributions and not dividing at the end yields inflated counts.\n• Wrong reflection types in D_n: For even n, reflections split into two conjugacy classes; using a single reflection pattern gives wrong answers. For odd n, each reflection has one fixed point and (n-1)/2 transpositions.\n• Confusing cycle length with number of cycles: In p_j^{a_j}, j is the cycle length and a_j is how many cycles of that length appear. Don’t swap them.\n• Mishandling weights: To count by the number of black items with a total of k colors, the correct substitution is p_j -> z^j + (k-1). Using p_j -> z + (k-1) ignores the fact that a j-cycle colored black contributes j black positions.\n• Overflow and non-integer division: Counts grow fast; use big integers. Intermediate sums before averaging must be exactly divisible by |G|; don’t truncate.\n• Double-counting reflections and rotations in D_n: The dihedral group has n rotations and n reflections. Mixing up 2n with n in averages causes off-by-2 errors.\n• Building the wrong action: Ensure your permutation truly describes how positions move under symmetry. For grids, define indices consistently and test that generators compose to the expected group size.

Key Formulas

Burnside's Lemma

Orbits = \frac{1}{∣ G ∣} g \in G \sum ∣ Fix (g) ∣

Explanation: The number of distinct colorings up to the group action equals the average number of colorings fixed by each group element. This is the foundation PET builds on.

Cycle Index Polynomial

Z (G) = \frac{1}{∣ G ∣} g \in G \sum j \geq 1 \prod p_{j}^{a_{j} (g)}

Explanation: Each group element contributes a monomial whose exponents are the counts of cycles of each length. Averaging over the group yields Z(G), a compact encoding of all cycle structures.

Pólya Substitution

F_{G} = Z (G)_{p_{j} \to \sum_{c \in C} w (c)^{j}} = \frac{1}{∣ G ∣} g \in G \sum j \geq 1 \prod (c \in C \sum w (c)^{j})^{a_{j} (g)}

Explanation: Replace each $p_{j}$ by the sum of j-th powers of color weights. The resulting polynomial (or number) counts weighted colorings up to symmetry.

Cycle Index of the Cyclic Group

Z (C_{n}) = \frac{1}{n} d ∣ n \sum φ (d) p_{d}^{n / d}

Explanation: Rotations of an n-gon decompose positions into cycles of equal length. Grouping by gcd and using Euler's totient gives this closed form.

Cycle Index of the Dihedral Group

Z (D_{n}) = {\frac{1}{2} Z (C_{n}) + \frac{1}{2} p_{1} p_{2}^{(n - 1) /2}, \frac{1}{2} Z (C_{n}) + \frac{1}{4} (p_{1}^{2} p_{2}^{n /2 - 1} + p_{2}^{n /2}), n odd, n even .

Explanation: Reflections have different cycle structures depending on parity of n. Averaging rotations and reflections yields these formulas.

Necklaces with k Colors

N_{necklace} (n, k) = \frac{1}{n} d ∣ n \sum φ (d) k^{n / d}

Explanation: Number of distinct colorings of an n-bead necklace with k colors up to rotation. It follows from substituting $p_{j}$ -> k into Z( $C_{n}$ ).

Bracelets with k Colors

B_{bracelet} (n, k) = Z (D_{n})_{p_{j} \to k}

Explanation: Number of bracelets (rotations and reflections) with k colors is obtained by evaluating the dihedral cycle index at $p_{j} = k$ .

Black-Count Enumerator

E_{G} (z) = Z (G)_{p_{j} \to z^{j} + (k - 1)}

Explanation: Tracks orbits by the number of black positions when there are k total colors and one distinguished color (black). The coefficient of $z^{b}$ equals the number of orbits with exactly b black positions.

Single-Length Cycle Factor

i = 0 \sum g (i g) (k - 1)^{g - i} z^{i ℓ}

Explanation: When an element has g cycles of equal length $ℓ$ , their contribution after substitution p_ $ℓ$ -> $z^{ℓ}$ + (k-1) is this binomial expansion.

Totient Summation Identity

d ∣ n \sum φ (d) = n

Explanation: Useful sanity check when summing over divisors in $C_{n}$ computations. The identity ensures correct averaging.

Complexity Analysis

Algorithmic complexity in PET comes from two sources: (1) enumerating/grouping group elements by cycle type, and (2) polynomial multiplications after substitution. For groups with closed-form cycle indices (

C_{n}

D_{n}

), enumeration is O(

τ

(n)) or O(n) depending on the formula used. For

C_{n}

necklaces with k colors, computing

N_{n} ec k l a ce

uses the divisor sum with Euler’s totient, which is O(

τ

(n) +

lo g

n) exponentiations; with fast exponentiation and factoring n to get divisors, this is typically O(

n

) preprocessing and O(

τ

(n)

lo g

n) arithmetic. Space is O(1) beyond storing divisors. For

D_{n}

bracelets with black-count distribution, rotations contribute polynomials with degrees up to n and structured gaps (only multiples of certain lengths). We can compute each rotation’s contribution via a binomial expansion in O(g) where

g = g c d (n

, r). Summing over

r = 0

..n-1 gives total O(

\sum_{r}

gcd(n,r)) = O(n

τ

(n)) arithmetic on big integers. Reflections add O(n) more operations with cheap binomial factors. The polynomial degree is n, so storing the accumulator uses O(n) big-integer coefficients. In general Burnside with arbitrary generators, generating the group of size

∣ G ∣

via closure on permutations of size n is O(

∣ G ∣

n) for composition and hashing, and then for each element we factor cycles in O(n). Multiplying cycle factors naively is O(

n^{2}

), but aggregating equal-length cycles and using binomial expansions reduces it to roughly O(n

n

) in typical small-n cases. Space is O(

∣ G ∣

n) to store permutations if retained; for on-the-fly enumeration (e.g., BFS with seen set) space is O(

∣ G ∣

). Big-integer arithmetic can dominate for larger n and k, so use efficient exponentiation and avoid dense convolutions when closed forms (binomial gaps) exist.

Code Examples

Necklaces under rotation (C_n): count with k colors using Euler’s totient

1 #include <bits/stdc++.h>
2 #include <boost/multiprecision/cpp_int.hpp>
3 using namespace std;
4 using boost::multiprecision::cpp_int;
5 
6 // Fast exponentiation for big integers
7 cpp_int ipow(cpp_int a, long long e) {
8     cpp_int r = 1;
9     while (e > 0) {
10         if (e & 1) r *= a;
11         a *= a;
12         e >>= 1;
13     }
14     return r;
15 }
16 
17 // Factor n and generate all divisors
18 vector<long long> divisors_of(long long n) {
19     vector<long long> small, large;
20     for (long long p = 1; p * p <= n; ++p) {
21         if (n % p == 0) {
22             small.push_back(p);
23             if (p * p != n) large.push_back(n / p);
24         }
25     }
26     reverse(large.begin(), large.end());
27     vector<long long> divs = small;
28     divs.insert(divs.end(), large.begin(), large.end());
29     return divs;
30 }
31 
32 // Compute phi(d) from its prime factors (which divide n)
33 long long phi_from_factors(long long d, const vector<long long>& primes) {
34     long long res = d;
35     for (long long p : primes) {
36         if (p * p > d) break;
37         if (d % p == 0) {
38             while (d % p == 0) d /= p;
39             res = res / p * (p - 1);
40         }
41     }
42     if (d > 1) res = res / d * (d - 1);
43     return res;
44 }
45 
46 // Prime factors of n (distinct)
47 vector<long long> distinct_prime_factors(long long n) {
48     vector<long long> primes;
49     for (long long p = 2; p * p <= n; ++p) {
50         if (n % p == 0) {
51             primes.push_back(p);
52             while (n % p == 0) n /= p;
53         }
54     }
55     if (n > 1) primes.push_back(n);
56     return primes;
57 }
58 
59 int main() {
60     ios::sync_with_stdio(false);
61     cin.tie(nullptr);
62 
63     long long n; long long k;
64     // Input: n beads, k colors
65     // Example: 12 3
66     if (!(cin >> n >> k)) return 0;
67 
68     vector<long long> divs = divisors_of(n);
69     vector<long long> primes = distinct_prime_factors(n);
70 
71     cpp_int total = 0;
72     for (long long d : divs) {
73         long long phi_d = phi_from_factors(d, primes);
74         // term: phi(d) * k^{n/d}
75         cpp_int term = cpp_int(phi_d) * ipow(cpp_int(k), n / d);
76         total += term;
77     }
78     // Answer: (1/n) * sum_{d|n} phi(d) * k^{n/d}
79     cpp_int ans = total / n; // exactly divisible
80 
81     cout << ans << "\n";
82     return 0;
83 }
84

Implements the closed-form necklace formula N(n,k) = (1/n) sum_{d|n} phi(d) k^{n/d}. We compute divisors of n, evaluate Euler’s totient for each divisor using primes of n, and use big-integer exponentiation to avoid overflow.

Time: O(√n + τ(n) log n) big-integer operationsSpace: O(τ(n)) for divisor storage

Bracelets (D_n) distribution by number of black beads among k colors

1 #include <bits/stdc++.h>
2 #include <boost/multiprecision/cpp_int.hpp>
3 using namespace std;
4 using boost::multiprecision::cpp_int;
5 
6 // Polynomial as coefficients of z^i (0..n)
7 using Poly = vector<cpp_int>;
8 
9 cpp_int ipow(cpp_int a, long long e){ cpp_int r=1; while(e){ if(e&1) r*=a; a*=a; e>>=1;} return r; }
10 
11 // Add polynomial b into a (resize as needed)
12 void poly_add(Poly &a, const Poly &b){
13     if (a.size() < b.size()) a.resize(b.size());
14     for (size_t i=0;i<b.size();++i) a[i]+=b[i];
15 }
16 
17 // Convolution (Cauchy) up to degree n
18 Poly poly_mul(const Poly &a, const Poly &b, int n){
19     Poly c(min((int)a.size()+(int)b.size()-1, n+1));
20     for (int i=0;i<(int)a.size();++i){
21         for (int j=0;j<(int)b.size();++j){
22             int d=i+j; if (d>n) break; c[d]+=a[i]*b[j];
23         }
24     }
25     return c;
26 }
27 
28 // Build (z^len + (k-1))^g as a sparse polynomial using binomial theorem
29 Poly binom_gap(long long len, long long g, long long km1, int n){
30     Poly p(n+1); // degree up to n
31     for (long long i=0;i<=g;i++){
32         // binomial(g, i) * (k-1)^{g-i} * z^{len*i}
33         // compute C(g,i) multiplicatively to avoid overflow
34     }
35     // We'll compute C(g,i) iteratively
36     cpp_int comb = 1; // C(g,0)
37     for (long long i=0;i<=g;i++){
38         cpp_int term = comb * ipow(cpp_int(km1), g - i);
39         long long deg = len * i;
40         if (deg <= n) p[deg] += term;
41         // update comb to C(g, i+1)
42         if (i < g) comb = comb * (g - i) / (i + 1);
43     }
44     return p;
45 }
46 
47 int main(){
48     ios::sync_with_stdio(false);
49     cin.tie(nullptr);
50 
51     long long n, k; // n beads, k colors total, track black among them
52     // Example: 10 3
53     if(!(cin >> n >> k)) return 0;
54     long long km1 = k - 1;
55 
56     // Accumulator for sum over all 2n group elements
57     Poly acc(n+1, 0);
58 
59     // Rotations r=0..n-1: cycle structure is g=gcd(n,r) cycles of length L=n/g
60     for (long long r=0; r<n; ++r){
61         long long g = std::gcd(n, r);
62         long long L = n / g;
63         Poly pr = binom_gap(L, g, km1, (int)n); // (z^L + (k-1))^g
64         poly_add(acc, pr);
65     }
66 
67     // Reflections: depends on parity of n
68     if (n % 2 == 1){
69         // n reflections, each: 1 fixed point (len 1), (n-1)/2 transpositions (len 2)
70         Poly p1(n+1); p1[0] = km1; if (n>=1) p1[1] += 1; // (z + (k-1))
71         Poly p2 = binom_gap(2, (n-1)/2, km1, (int)n); // (z^2 + (k-1))^{(n-1)/2}
72         Poly pref = poly_mul(p1, p2, (int)n);
73         for (int t=0;t<(int)n;t++){} // no-op to stress bounds
74         // add n copies
75         for (long long i=0;i<n;i++) poly_add(acc, pref);
76     } else {
77         // even n: n/2 reflections of type A: 2 fixed points + (n/2-1) transpositions
78         // (z + (k-1))^2 * (z^2 + (k-1))^{n/2 - 1}
79         Poly p1( min((int)n+1, 3) ); // (z + (k-1))^2 = z^2 + 2(k-1)z + (k-1)^2
80         if ((int)p1.size() < 3) p1.resize(3);
81         p1[0] = cpp_int(km1) * cpp_int(km1);
82         p1[1] = cpp_int(2) * cpp_int(km1);
83         p1[2] = 1;
84         Poly p2 = binom_gap(2, n/2 - 1, km1, (int)n);
85         Poly prefA = poly_mul(p1, p2, (int)n);
86         for (long long i=0;i<n/2;i++) poly_add(acc, prefA);
87         // n/2 reflections of type B: all transpositions (no fixed points): (z^2 + (k-1))^{n/2}
88         Poly prefB = binom_gap(2, n/2, km1, (int)n);
89         for (long long i=0;i<n/2;i++) poly_add(acc, prefB);
90     }
91 
92     // Average over |D_n| = 2n
93     for (int d=0; d<=n; ++d) acc[d] /= (2*n);
94 
95     // Output: counts for exactly b black beads, b=0..n
96     for (int b=0; b<=n; ++b){
97         if (b) cout << ' ';
98         cout << acc[b];
99     }
100     cout << "\n";
101     return 0;
102 }
103

Computes the z-enumerator E(z) = (1/|D_n|) sum_g ∏_cycles (z^{len} + (k-1)), whose z^b coefficient is the number of bracelets with exactly b black beads (and k−1 other colors). Rotations are handled via binomial expansions with gaps. Reflections split by parity of n, using the correct cycle structures. We average over 2n group elements and print coefficients 0..n.

Time: O(n τ(n)) big-integer ops for rotations plus O(n) for reflections; total polynomial degree nSpace: O(n) big-integer coefficients

General PET via generators: 2×2 grid under D4 (rotations/reflections), distribution of black cells among k colors

1 #include <bits/stdc++.h>
2 #include <boost/multiprecision/cpp_int.hpp>
3 using namespace std;
4 using boost::multiprecision::cpp_int;
5 
6 using Perm = vector<int>; // permutation of {0..n-1}
7 using Poly = vector<cpp_int>; // coefficients in z
8 
9 cpp_int ipow(cpp_int a, long long e){ cpp_int r=1; while(e){ if(e&1) r*=a; a*=a; e>>=1;} return r; }
10 
11 Perm compose(const Perm &a, const Perm &b){ // a∘b (apply b then a)
12     int n = (int)a.size();
13     Perm c(n);
14     for (int i=0;i<n;++i) c[i] = a[b[i]];
15     return c;
16 }
17 
18 struct PermHash {
19     size_t operator()(const Perm& p) const noexcept {
20         size_t h = 1469598103934665603ull;
21         for (int x: p){ h ^= (size_t)x + 0x9e3779b97f4a7c15ull + (h<<6) + (h>>2); }
22         return h;
23     }
24 };
25 
26 // Find cycles of permutation p on {0..n-1}
27 vector<int> cycle_lengths(const Perm &p){
28     int n = (int)p.size();
29     vector<int> vis(n,0), lens;
30     for (int i=0;i<n;++i){
31         if (!vis[i]){
32             int v=i, len=0; while(!vis[v]){ vis[v]=1; v=p[v]; ++len; }
33             lens.push_back(len);
34         }
35     }
36     return lens;
37 }
38 
39 // Multiply polynomials up to degree n
40 Poly poly_mul(const Poly &a, const Poly &b, int n){
41     Poly c(min((int)a.size()+(int)b.size()-1, n+1));
42     for (int i=0;i<(int)a.size();++i){
43         for (int j=0;j<(int)b.size();++j){
44             int d=i+j; if (d>n) break; c[d]+=a[i]*b[j];
45         }
46     }
47     return c;
48 }
49 
50 // Build (z^len + (k-1))^t quickly (binomial with gaps)
51 Poly binom_gap(int len, int t, int km1, int n){
52     Poly p(n+1);
53     cpp_int comb=1; // C(t,0)
54     for (int i=0;i<=t;++i){
55         cpp_int term = comb * ipow(cpp_int(km1), t - i);
56         int deg = len * i; if (deg <= n) p[deg] += term;
57         if (i < t) comb = comb * (t - i) / (i + 1);
58     }
59     return p;
60 }
61 
62 int main(){
63     ios::sync_with_stdio(false);
64     cin.tie(nullptr);
65 
66     int k; // total colors; track black among them
67     // Example input: 3
68     if(!(cin >> k)) return 0;
69     int km1 = k - 1;
70 
71     // 2x2 grid cells indexed as (r,c) -> r*2 + c
72     auto idx = [](int r,int c){ return r*2 + c; };
73 
74     // Generator g1: rotation by 90° clockwise: (r,c) -> (c, 1-r)
75     Perm rot4(4); for (int r=0;r<2;++r) for (int c=0;c<2;++c){ int from=idx(r,c); int nr=c, nc=1-r; rot4[from]=idx(nr,nc);}    
76     // Generator g2: reflection across vertical axis: (r,c) -> (r, 1-c)
77     Perm refl(4); for (int r=0;r<2;++r) for (int c=0;c<2;++c){ int from=idx(r,c); int nr=r, nc=1-c; refl[from]=idx(nr,nc);}   
78 
79     // Generate the group <rot4, refl>
80     vector<Perm> gens = {rot4, refl};
81     unordered_set<Perm, PermHash> seen;
82     queue<Perm> q;
83 
84     // identity
85     Perm id(4); iota(id.begin(), id.end(), 0);
86     seen.insert(id); q.push(id);
87     while(!q.empty()){
88         Perm cur = q.front(); q.pop();
89         for (const Perm &g: gens){
90             Perm nxt = compose(g, cur);
91             if (!seen.count(nxt)){ seen.insert(nxt); q.push(nxt); }
92         }
93     }
94     vector<Perm> group(seen.begin(), seen.end());
95     // Expect |D4| = 8
96 
97     int n = 4; // positions
98     Poly acc(n+1); // sum of polynomials over group
99     for (const Perm &p : group){
100         // Count multiplicities of cycle lengths
101         vector<int> lens = cycle_lengths(p);
102         // Build factor per length: (z^L + (k-1))^{count}
103         // Aggregate equal lengths to reduce multiplications
104         unordered_map<int,int> cnt;
105         for(int L: lens) cnt[L]++;
106         Poly poly(1, cpp_int(1)); // start with 1
107         for (auto [L, t] : cnt){
108             Poly factor = binom_gap(L, t, km1, n);
109             poly = poly_mul(poly, factor, n);
110         }
111         // add to accumulator
112         if (acc.size() < poly.size()) acc.resize(poly.size());
113         for (size_t d=0; d<poly.size(); ++d) acc[d] += poly[d];
114     }
115 
116     // Average
117     cpp_int Gsize = (cpp_int)group.size();
118     for (int d=0; d<=n; ++d) acc[d] /= Gsize;
119 
120     // Output distribution: number of inequivalent colorings by exactly b black cells
121     for (int b=0;b<=n;++b){ if (b) cout << ' '; cout << acc[b]; }
122     cout << "\n";
123 
124     return 0;
125 }
126

Generates the full D4 symmetry group on the 2×2 grid from two intuitive generators (90° rotation and vertical reflection), then applies Burnside/PET with the substitution p_j -> z^j + (k−1). Cycle decompositions are computed for each permutation, factors are multiplied using binomial expansions, and coefficients give the orbit counts by black-cell count.

Time: O(|G| · n + |G| · P(n)) where P(n) is polynomial multiplication cost up to degree n (here n=4, trivial). Group generation is O(|G|·n).Space: O(|G|·n) to store permutations; O(n) for the accumulator polynomial

1	#include <bits/stdc++.h>
2	#include <boost/multiprecision/cpp_int.hpp>
3	using namespace std;
4	using boost::multiprecision::cpp_int;
5
6	// Fast exponentiation for big integers
7	cpp_int ipow(cpp_int a, long long e) {
8	cpp_int r = 1;
9	while (e > 0) {
10	if (e & 1) r *= a;
11	a *= a;
12	e >>= 1;
13	}
14	return r;
15	}
16
17	// Factor n and generate all divisors
18	vector<long long> divisors_of(long long n) {
19	vector<long long> small, large;
20	for (long long p = 1; p * p <= n; ++p) {
21	if (n % p == 0) {
22	small.push_back(p);
23	if (p * p != n) large.push_back(n / p);
24	}
25	}
26	reverse(large.begin(), large.end());
27	vector<long long> divs = small;
28	divs.insert(divs.end(), large.begin(), large.end());
29	return divs;
30	}
31
32	// Compute phi(d) from its prime factors (which divide n)
33	long long phi_from_factors(long long d, const vector<long long>& primes) {
34	long long res = d;
35	for (long long p : primes) {
36	if (p * p > d) break;
37	if (d % p == 0) {
38	while (d % p == 0) d /= p;
39	res = res / p * (p - 1);
40	}
41	}
42	if (d > 1) res = res / d * (d - 1);
43	return res;
44	}
45
46	// Prime factors of n (distinct)
47	vector<long long> distinct_prime_factors(long long n) {
48	vector<long long> primes;
49	for (long long p = 2; p * p <= n; ++p) {
50	if (n % p == 0) {
51	primes.push_back(p);
52	while (n % p == 0) n /= p;
53	}
54	}
55	if (n > 1) primes.push_back(n);
56	return primes;
57	}
58
59	int main() {
60	ios::sync_with_stdio(false);
61	cin.tie(nullptr);
62
63	long long n; long long k;
64	// Input: n beads, k colors
65	// Example: 12 3
66	if (!(cin >> n >> k)) return 0;
67
68	vector<long long> divs = divisors_of(n);
69	vector<long long> primes = distinct_prime_factors(n);
70
71	cpp_int total = 0;
72	for (long long d : divs) {
73	long long phi_d = phi_from_factors(d, primes);
74	// term: phi(d) * k^{n/d}
75	cpp_int term = cpp_int(phi_d) * ipow(cpp_int(k), n / d);
76	total += term;
77	}
78	// Answer: (1/n) * sum_{d\|n} phi(d) * k^{n/d}
79	cpp_int ans = total / n; // exactly divisible
80
81	cout << ans << "\n";
82	return 0;
83	}
84