Enumerative combinatorics

"""Integer partitions: generating functions, Euler's theorem and Glaisher's map.

We expand the generating functions

    prod 1/(1 - x^k)          (all partitions, p(n)),
    prod 1/(1 - x^(2k-1))     (partitions into odd parts),
    prod (1 + x^k)            (partitions into distinct parts)

as exact integer power series up to x^60, and compare the last two, as
Euler's theorem predicts. For n = 8 we also list Glaisher's bijection, which
merges equal odd parts into distinct parts, and check that it is a bijection.
Finally we compare p(n) with the Hardy-Ramanujan asymptotic formula and test
Ramanujan's congruence p(5m + 4) = 0 (mod 5) on a range.

All counts are exact. Agreement up to x^60 is evidence; Euler's theorem is
proved for all n in the lesson (and in Mathlib).
"""

import math

def times_geometric(series, step):
    """Multiply a truncated power series by 1/(1 - x^step), in place.

    If b = a / (1 - x^step) then b = a + x^step * b, so b[i] = a[i] + b[i - step].
    """
    for i in range(step, len(series)):
        series[i] += series[i - step]
    return series

def times_one_plus(series, step):
    """Multiply a truncated power series by (1 + x^step), in place."""
    for i in range(len(series) - 1, step - 1, -1):
        series[i] += series[i - step]
    return series

def partition_counts(n_max):
    """p(0), ..., p(n_max): coefficients of the product of 1/(1 - x^k)."""
    series = [1] + [0] * n_max
    for k in range(1, n_max + 1):
        times_geometric(series, k)
    return series

def odd_part_counts(n_max):
    """Coefficients of the product of 1/(1 - x^k) over odd k."""
    series = [1] + [0] * n_max
    for k in range(1, n_max + 1, 2):
        times_geometric(series, k)
    return series

def distinct_part_counts(n_max):
    """Coefficients of the product of (1 + x^k)."""
    series = [1] + [0] * n_max
    for k in range(1, n_max + 1):
        times_one_plus(series, k)
    return series

def partitions(n, largest=None):
    """All partitions of n as non-increasing tuples."""
    if largest is None:
        largest = n
    if n == 0:
        return [()]
    out = []
    for first in range(min(n, largest), 0, -1):
        out += [(first,) + rest for rest in partitions(n - first, first)]
    return out

def glaisher(odd_partition):
    """Merge equal odd parts: m copies of k become k * 2^i for the bits of m."""
    parts = []
    for k in set(odd_partition):
        m = odd_partition.count(k)
        i = 0
        while m:
            if m & 1:
                parts.append(k * 2 ** i)
            m >>= 1
            i += 1
    return tuple(sorted(parts, reverse=True))

def hardy_ramanujan(n):
    """The leading asymptotic exp(pi sqrt(2n/3)) / (4 n sqrt 3)."""
    return math.exp(math.pi * math.sqrt(2 * n / 3)) / (4 * n * math.sqrt(3))

def main():
    n_max = 60
    odd = odd_part_counts(n_max)
    distinct = distinct_part_counts(n_max)
    print("Generating functions expanded exactly up to x^60:")
    print(f"  odd-part counts:      {odd[:13]} ...")
    print(f"  distinct-part counts: {distinct[:13]} ...")
    print(f"  equal for every n <= {n_max}: {odd == distinct}")
    print()
    eights = partitions(8)
    odds8 = [q for q in eights if all(part % 2 for part in q)]
    images = [glaisher(q) for q in odds8]
    print(f"Glaisher's map on the {len(odds8)} odd-part partitions of 8:")
    for q, image in zip(odds8, images):
        print(f"  {'+'.join(map(str, q)):<17} -> {'+'.join(map(str, image))}")
    distinct8 = {q for q in eights if len(set(q)) == len(q)}
    print(f"  images distinct and exactly the distinct-part partitions: "
          f"{len(set(images)) == len(images) and set(images) == distinct8}")
    print()
    p = partition_counts(300)
    print("p(n) versus Hardy-Ramanujan exp(pi sqrt(2n/3))/(4n sqrt 3):")
    for n in (10, 50, 100, 300):
        print(f"  n = {n:>3}: p(n) = {p[n]:<19} ratio = "
              f"{p[n] / hardy_ramanujan(n):.4f}")
    bad = [n for n in range(4, 300, 5) if p[n] % 5]
    print(f"Ramanujan: p(5m + 4) divisible by 5 for all 5m + 4 < 300: {not bad}")
    print()
    print("Epistemic status: exact integer computation, exhaustive for n <= 60 "
          "(Euler)\nand n < 300 (Ramanujan), and for n = 8 (Glaisher). These "
          "are finite checks\nof theorems proved in general; the asymptotic "
          "ratios are numerical and only\nillustrate that the relative error "
          "tends to 0, slowly.")

if __name__ == "__main__":
    main()

"""Exponential generating functions, computed exactly and checked by brute force.

A labelled structure on n points is counted by a_n, and its exponential
generating function (EGF) is A(x) = sum a_n x^n / n!. The exponential
formula says: if a structure is a set of "connected" pieces counted by
C(x), then A(x) = exp(C(x)). We expand three EGFs with exact fractions:

* set partitions:   exp(e^x - 1)        (blocks are nonempty sets)
* involutions:      exp(x + x^2/2)      (cycles of length 1 or 2)
* derangements:     e^(-x) / (1 - x)    (permutations with no fixed point)

and compare their coefficients with direct enumeration for n <= 7. We also
watch D(n)/n! approach 1/e, with the alternating-series error bound.
"""

import itertools
import math
from fractions import Fraction

N = 12

def series_mul(a, b):
    """Product of two truncated power series (lists of Fractions)."""
    return [sum(a[i] * b[n - i] for i in range(n + 1)) for n in range(len(a))]

def series_exp(f):
    """exp(f) for a truncated series with f[0] = 0, via g' = f' g:
    n g_n = sum_{k=1}^{n} k f_k g_{n-k}."""
    if f[0] != 0:
        raise ValueError("constant term must be 0")
    g = [Fraction(1)] + [Fraction(0)] * (len(f) - 1)
    for n in range(1, len(f)):
        g[n] = sum(k * f[k] * g[n - k] for k in range(1, n + 1)) / n
    return g

def exp_series(c=1):
    """e^(c x) as a truncated series."""
    return [Fraction(c) ** n / math.factorial(n) for n in range(N + 1)]

def egf_counts(series):
    """Turn EGF coefficients into counts a_n = n! [x^n]."""
    counts = [s * math.factorial(n) for n, s in enumerate(series)]
    assert all(c.denominator == 1 for c in counts)
    return [int(c) for c in counts]

def set_partitions_egf():
    e = exp_series()
    return egf_counts(series_exp([Fraction(0)] + e[1:]))

def involutions_egf():
    f = [Fraction(0)] * (N + 1)
    f[1], f[2] = Fraction(1), Fraction(1, 2)
    return egf_counts(series_exp(f))

def derangements_egf():
    geometric = [Fraction(1)] * (N + 1)
    return egf_counts(series_mul(exp_series(-1), geometric))

def count_set_partitions(n):
    """Brute force: build partitions of {0..n-1} element by element."""
    blocks_list = [[]]
    for x in range(n):
        new = []
        for blocks in blocks_list:
            for i in range(len(blocks)):
                new.append(blocks[:i] + [blocks[i] + [x]] + blocks[i + 1:])
            new.append(blocks + [[x]])
        blocks_list = new
    return len(blocks_list)

def count_involutions(n):
    return sum(1 for p in itertools.permutations(range(n))
               if all(p[p[i]] == i for i in range(n)))

def count_derangements(n):
    return sum(1 for p in itertools.permutations(range(n))
               if all(p[i] != i for i in range(n)))

def main():
    rows = [
        ("set partitions (Bell)", "exp(e^x - 1)", set_partitions_egf(),
         count_set_partitions),
        ("involutions", "exp(x + x^2/2)", involutions_egf(), count_involutions),
        ("derangements", "e^(-x)/(1 - x)", derangements_egf(),
         count_derangements),
    ]
    for name, egf, coeffs, brute in rows:
        checked = all(coeffs[n] == brute(n) for n in range(8))
        print(f"{name} from {egf}:")
        print(f"  n = 0..10: {coeffs[:11]}")
        print(f"  equal to brute-force counts for n <= 7: {checked}")
    print()
    d = derangements_egf()
    print("Derangements: D(n)/n! approaches 1/e = 0.3678794412...")
    for n in (3, 6, 9, 12):
        ratio = Fraction(d[n], math.factorial(n))
        error = abs(float(ratio) - math.exp(-1))
        bound = 1 / math.factorial(n + 1)
        print(f"  n = {n:>2}: D(n)/n! = {float(ratio):.10f}, error {error:.2e}"
              f" <= {bound:.2e}: {error <= bound}")
    print(f"  D(n) = round(n!/e) for 1 <= n <= 12: "
          f"{all(d[n] == round(math.factorial(n) / math.e) for n in range(1, 13))}")
    print()
    print("Epistemic status: exact rational arithmetic for the series and "
          "exhaustive\nenumeration for n <= 7, so the agreement is a finite "
          "check of the\nexponential formula, not a proof; the error column "
          "uses floating point,\nwhich is ample at this precision.")

if __name__ == "__main__":
    main()

"""Counting up to symmetry: necklaces, bracelets and coloured cubes.

Two colourings are "the same" if a symmetry carries one to the other, so we
are counting orbits of a group acting on colourings. We count orbits in two
independent ways and compare:

* brute force: list all colourings and pick one canonical representative
  (the smallest) from each orbit;
* Burnside's lemma: average the number of colourings fixed by each group
  element, which for necklaces gives (1/n) sum_{d | n} phi(d) k^(n/d), and
  for the cube is Polya's cycle-index polynomial evaluated at k colours.

We also show why the naive count k^n / n is wrong: symmetric colourings
have smaller orbits. The checks are exhaustive for the sizes shown.
"""

import itertools
from fractions import Fraction
from math import gcd

def phi(n):
    """Euler's totient function."""
    return sum(1 for a in range(1, n + 1) if gcd(a, n) == 1)

def rotations(n):
    """The n rotations of an n-bead necklace, as position permutations."""
    return [tuple((i + r) % n for i in range(n)) for r in range(n)]

def dihedral(n):
    """Rotations and reflections: the dihedral group of order 2n."""
    return rotations(n) + [tuple((r - i) % n for i in range(n)) for r in range(n)]

def act(perm, colouring):
    """The colouring moved by a permutation of positions."""
    return tuple(colouring[perm[i]] for i in range(len(perm)))

def orbits_brute_force(group, k, size):
    """Number of orbits: count colourings that are smallest in their orbit."""
    return sum(1 for c in itertools.product(range(k), repeat=size)
               if all(c <= act(g, c) for g in group))

def cycles(perm):
    """Number of cycles of a permutation."""
    seen, count = set(), 0
    for start in range(len(perm)):
        if start not in seen:
            count += 1
            i = start
            while i not in seen:
                seen.add(i)
                i = perm[i]
    return count

def orbits_burnside(group, k):
    """Burnside: average of |Fix g|; a permutation with c cycles fixes k^c."""
    total = sum(k ** cycles(g) for g in group)
    return Fraction(total, len(group))

def necklace_formula(n, k):
    """(1/n) sum over divisors d of n of phi(d) k^(n/d)."""
    return Fraction(sum(phi(d) * k ** (n // d) for d in range(1, n + 1)
                        if n % d == 0), n)

def cube_rotations():
    """The 24 rotations of a cube acting on its faces U, D, F, B, L, R."""
    spin = (0, 1, 5, 4, 2, 3)   # quarter turn about the U-D axis
    tilt = (2, 3, 1, 0, 4, 5)   # quarter turn about the L-R axis
    group, frontier = {tuple(range(6))}, [tuple(range(6))]
    while frontier:
        g = frontier.pop()
        for h in (spin, tilt):
            gh = tuple(g[h[i]] for i in range(6))
            if gh not in group:
                group.add(gh)
                frontier.append(gh)
    return sorted(group)

def main():
    print("Necklaces of n beads in k colours (rotations only):")
    print("   n  k   brute force   Burnside   formula   naive k^n/n")
    for n, k in [(4, 2), (6, 2), (6, 3), (8, 2), (9, 2), (9, 3)]:
        group = rotations(n)
        brute = orbits_brute_force(group, k, n)
        burn = orbits_burnside(group, k)
        naive = Fraction(k ** n, n)
        print(f"{n:>4} {k:>2} {brute:>13} {str(burn):>10} "
              f"{str(necklace_formula(n, k)):>9} {float(naive):>13.2f}")
    print("Bracelets (rotations and reflections), k = 2:")
    for n in (4, 6, 8, 9):
        group = dihedral(n)
        print(f"   n = {n}: brute force {orbits_brute_force(group, 2, n)}, "
              f"Burnside {orbits_burnside(group, 2)}")
    print()
    cube = cube_rotations()
    types = sorted(cycles(g) for g in cube)
    print(f"Cube: {len(cube)} rotations generated from two quarter turns.")
    print(f"  number with c cycles on the faces, c = 1..6: "
          f"{[types.count(c) for c in range(1, 7)]}")
    for k in (2, 3):
        print(f"  face colourings with {k} colours up to rotation: brute force "
              f"{orbits_brute_force(cube, k, 6)}, Burnside "
              f"{orbits_burnside(cube, k)}")
    print()
    print("Epistemic status: exact and exhaustive for every size shown - each "
          "count\nwas obtained by listing all colourings and agrees with "
          "Burnside's lemma.\nThat is evidence for the formulas at these sizes; "
          "the lemma itself is\nproved in the lesson (and in Mathlib).")

if __name__ == "__main__":
    main()