Necklaces, bracelets and cubes: Burnside versus brute force

"""Counting patterns up to symmetry: Burnside's lemma against brute force.

A necklace is a ring of n coloured beads up to rotation; a bracelet is the same up
to rotation and reflection (the dihedral group). Burnside's lemma counts orbits as
the average number of fixed colourings, and a permutation with c cycles fixes
exactly k**c colourings with k colours. The program compares that formula with a
brute-force count of orbits (one canonical representative per orbit) for small n,
and does the same for colourings of the six faces of a cube under its 24
rotations, which it generates from two quarter turns. Exact and exhaustive for
the sizes shown.
"""

from itertools import product

def compose(p, q):
    """Return the permutation 'first q, then p'."""
    return tuple(p[i] for i in q)

def closure(gens):
    """The permutation group generated by gens (all of the same degree)."""
    identity = tuple(range(len(gens[0])))
    elements, frontier = {identity}, [identity]
    while frontier:
        new = []
        for x in frontier:
            for g in gens:
                y = compose(g, x)
                if y not in elements:
                    elements.add(y)
                    new.append(y)
        frontier = new
    return elements

def cyclic_group(n):
    """Rotations of n beads: i -> i + k (mod n)."""
    return [tuple((i + k) % n for i in range(n)) for k in range(n)]

def dihedral_group(n):
    """Rotations and reflections of n beads: i -> +-i + k (mod n)."""
    reflections = [tuple((k - i) % n for i in range(n)) for k in range(n)]
    return sorted(set(cyclic_group(n) + reflections))

def cycle_count(p):
    """Number of cycles of the permutation p (fixed points count as cycles)."""
    seen, cycles = set(), 0
    for start in range(len(p)):
        if start not in seen:
            cycles += 1
            i = start
            while i not in seen:
                seen.add(i)
                i = p[i]
    return cycles

def burnside(group, k):
    """Number of orbits on k-colourings: average of k**cycles over the group."""
    total = sum(k ** cycle_count(g) for g in group)
    if total % len(group):
        raise ArithmeticError("Burnside sum not divisible by |G|")
    return total // len(group)

def brute_force(group, k):
    """Count orbits directly: keep the least colouring in each orbit."""
    n = len(group[0])
    reps = set()
    for colouring in product(range(k), repeat=n):
        images = (tuple(colouring[g[i]] for i in range(n)) for g in group)
        reps.add(min(images))
    return len(reps)

def cube_rotations():
    """The 24 rotations of a cube as permutations of its faces U D F B L R."""
    quarter_vertical = (0, 1, 5, 4, 2, 3)   # F->R->B->L->F, U and D fixed
    quarter_front = (5, 4, 2, 3, 0, 1)      # U->R->D->L->U, F and B fixed
    return sorted(closure([quarter_vertical, quarter_front]))

def main():
    print("Orbits = average over the group of k^(number of cycles)")
    print("beads  necklaces k=2  bracelets k=2  necklaces k=3  bracelets k=3")
    agree = True
    for n in range(1, 9):
        row = []
        for k in (2, 3):
            for group in (cyclic_group(n), dihedral_group(n)):
                value = burnside(group, k)
                agree &= value == brute_force(group, k)
                row.append(value)
        print(f"{n:>5}  {row[0]:>13}  {row[1]:>13}  {row[2]:>13}  {row[3]:>13}")
    print(f"Burnside agrees with brute force in every case checked: {agree}")
    cube = cube_rotations()
    print(f"Cube: {len(cube)} rotations generated from two quarter turns")
    for k in (2, 3):
        print(f"  face colourings with {k} colours up to rotation: "
              f"Burnside {burnside(cube, k)}, brute force {brute_force(cube, k)}")
    print()
    print("Epistemic status: exact, exhaustive counts for the sizes shown.")
    print("Agreement here is evidence; the proof that Burnside's formula")
    print("always counts orbits is in the lesson, by counting the pairs (g, x)")
    print("with g . x = x in two ways.")

if __name__ == "__main__":
    main()