Pólya counting

module Polya

Pólya / Burnside math primitives used by Enumlib's Enumlib.count_inequivalent. Independently usable for general group-action orbit counting.

Exports:

• polya_count(perm_group, k_or_mults) — Burnside-averaged orbit count.
• cycle_structure(permutation) — cycle lengths of a single permutation.
• aperiodic_orbit_count(perm_group, T_factors, k_or_mults) — Möbius-corrected primitive (aperiodic) orbit count.

See research.md §4.6 (Rosenbrock 2016) for the algorithm digest and research.md §5.2.1 for the super-periodicity policy that motivates the aperiodic variant.

polya_count(perm_group, k::Integer) -> BigInt

Burnside-averaged count of orbits of k-colorings under perm_group — the raw orbit count, including super-periodic orbits. For each ρ ∈ perm_group with c(ρ) cycles, fix(ρ) = k^c(ρ); orbit count = (1/|G|) Σ_ρ k^c(ρ).

Returns BigInt. Cost: O(|G| · n).

The aperiodic version is aperiodic_orbit_count.

Examples

Cyclic C_4 on 4 positions (the rotations of a 1×1×4 supercell), binary colorings (k = 2): Burnside sum is (2^4 + 2^1 + 2^2 + 2^1) / 4 = 24 / 4 = 6.

julia> c4 = [[1,2,3,4], [2,3,4,1], [3,4,1,2], [4,1,2,3]];

julia> polya_count(c4, 2)
6

The 6 orbits are: {AAAA}, {BBBB}, {AAAB-rotation}, {ABBB-rotation}, {AABB-rotation} (period 4), and {ABAB, BABA} (period 2).

polya_count(perm_group, multiplicities::AbstractVector{<:Integer}) -> BigInt

Burnside-averaged orbit count at fixed multiplicities. For each ρ with cycle lengths [c_1, ..., c_t], the number of fixed colorings is the coefficient of x_1^{a_1}⋯x_k^{a_k} in ∏_j (x_1^{c_j} + ⋯ + x_k^{c_j}) (HF 2012 §A.2). Computed via dynamic programming over partial-multiplicity states.

Cost per permutation: O(t · ∏(ai + 1) · k). Across the whole group: O(|G| · t · ∏(ai + 1) · k).

Examples

Cyclic C_4 with multiplicities [2, 2] (2 of each species in 4 positions). Hand-verifiable: 6 raw colorings (C(4, 2)) group into 2 orbits — {AABB, ABBA, BBAA, BAAB} (period 4) and {ABAB, BABA} (period 2).

julia> c4 = [[1,2,3,4], [2,3,4,1], [3,4,1,2], [4,1,2,3]];

julia> polya_count(c4, [2, 2])
2

polya_count(perm_group, allowed_labels::AbstractVector{BitSet}) -> BigInt
polya_count(perm_group, allowed_labels::AbstractVector{BitSet}, multiplicities) -> BigInt

Burnside-averaged orbit count for label-restricted positions. allowed_labels[p] is the set of labels (0-based, as in Enumlib.Site) permitted at supercell position p.

A coloring fixed by ρ is constant on each cycle of ρ, so a cycle can only carry a label allowed at every position it visits: fix(ρ) = ∏_cycles |⋂ allowed_labels|, and 0 when some cycle's intersection is empty. With a uniform allowed_labels this reproduces the scalar-k method exactly.

This is the count that matches length(enumerate(...)) whenever Sites are heterogeneous — zinc-blende, half/full-Heusler, perovskite, or any site pinned to one species. The scalar-k method assumes every position may take any of the k labels, and overcounts by orders of magnitude when that is false.

cycle_structure(perm) -> Vector{Int}

Cycle-length multiset of a permutation, sorted descending.

Examples

julia> cycle_structure([2, 1, 4, 3, 5])    # cycles (1 2)(3 4)(5)
3-element Vector{Int64}:
 2
 2
 1

julia> cycle_structure([2, 3, 4, 1])       # cycles (1 2 3 4) — one 4-cycle
1-element Vector{Int64}:
 4

aperiodic_orbit_count(perm_group, snf_diagonal, k_or_mults) -> BigInt

Aperiodic orbit count — orbits whose stabilizer in the supercell's translation subgroup T is trivial. Matches length(enumerate(...; include_superperiodic = false)) per supercell.

snf_diagonal is a 3-tuple (d_1, d_2, d_3) from the supercell HNF's Smith Normal Form, with d_1 | d_2 | d_3 (so T = Z/d_1 × Z/d_2 × Z/d_3 and |T| = d_1·d_2·d_3 = n). k_or_mults is either an Integer (unrestricted, k colors) or AbstractVector{<:Integer} (fixed multiplicities).

Burnside on the aperiodic-labeling subspace, with the per-permutation count expanded by Möbius inversion over the subgroup lattice of T:

N_aperiodic = (1/|G|) · Σ_{g ∈ G} Σ_{H ⊆ T} μ_T(H) · |fix(g) ∩ fix(H)|

where |fix(g) ∩ fix(H)| is the number of labelings fixed by both g and every element of H — equivalently, labelings constant on each joint orbit of {g} ∪ H_perms acting on supercell positions. For unrestricted k-colorings this is k^(num joint orbits); for fixed multiplicities it's the polynomial-coefficient dynamic program on the joint-orbit partition.

This formulation is correct for non-normal H (subgroups of T need not be normal in G — e.g., for SNF (1, 2, 2) where T = Z/2 × Z/2 is non-cyclic). Earlier sketches assumed H-normality and gave wrong counts on non-cyclic T.

At typical supercell sizes (|T| ≤ ~32), subgroup enumeration is cheap; the inner per-(g, H) joint-orbit computation is O(n + |H|) via union-find.

Examples

Same cyclic C_4 case as polya_count, showing the aperiodic-vs-raw contrast. The SNF diagonal (1, 1, 4) says T = Z/4 (the 1×1×4 supercell's translation subgroup).

Unrestricted binary (k=2): polya_count gives 6 orbits total; aperiodic count drops the 1-cycle orbits {AAAA}, {BBBB} and the period-2 orbit {ABAB, BABA}, leaving 3 period-4 orbits.

julia> c4 = [[1,2,3,4], [2,3,4,1], [3,4,1,2], [4,1,2,3]];

julia> aperiodic_orbit_count(c4, (1, 1, 4), 2)
3

Fixed-multiplicity [2, 2]: polya_count gives 2 orbits; aperiodic drops the period-2 {ABAB, BABA}, leaving 1 period-4 orbit {AABB, ABBA, BBAA, BAAB}.

julia> c4 = [[1,2,3,4], [2,3,4,1], [3,4,1,2], [4,1,2,3]];

julia> aperiodic_orbit_count(c4, (1, 1, 4), [2, 2])
1

aperiodic_orbit_count(perm_group, snf_diagonal, allowed_labels::AbstractVector{BitSet}) -> BigInt
aperiodic_orbit_count(perm_group, snf_diagonal, allowed_labels::AbstractVector{BitSet}, multiplicities) -> BigInt

Label-restricted counterpart of aperiodic_orbit_count. Identical Möbius inversion over the subgroup lattice of T; the per-(g, H) fixed-count intersects allowed_labels across each joint orbit instead of assuming all k labels everywhere.

Use these whenever the Sites are heterogeneous — see the restricted polya_count.