counting.py

1import sys2from fractions import Fraction3from math import factorial, gcd4from collections import Counter56sys.set_int_max_str_digits(1000000)78# ANALYTIC COUNTING - THE UNIVERSE WITHOUT ENUMERATION910def _mobius(n):11    x = n12    mu = 113    p = 214    while p * p <= x:15        if x % p == 0:16            x //= p17            if x % p == 0:18                return 019            mu = -mu20        p += 121    if x > 1:22        mu = -mu23    return mu2425def _divisors(n):26    return [d for d in range(1, n + 1) if n % d == 0]2728def _pos_block_cycles(length):29    out = Counter()30    for period in _divisors(length):31        strings = sum(_mobius(period // d) * 2 ** d for d in _divisors(period))32        out[period] += strings // period33    return out3435def _neg_block_cycles(length):36    def step(x):37        bits = [(x >> i) & 1 for i in range(length)]38        new = [bits[(i - 1) % length] for i in range(length)]39        new[0] ^= 140        return sum(new[i] << i for i in range(length))41    seen = [False] * (2 ** length)42    out = Counter()43    for start in range(2 ** length):44        if seen[start]:45            continue46        run = 047        j = start48        while not seen[j]:49            seen[j] = True50            j = step(j)51            run += 152        out[run] += 153    return out5455def _combine(c1, c2):56    out = Counter()57    for l1, n1 in c1.items():58        for l2, n2 in c2.items():59            g = gcd(l1, l2)60            out[l1 * l2 // g] += n1 * n2 * g61    return out6263def _class_cycles(pos, neg):64    blocks = [_pos_block_cycles(L) for L in pos] + [_neg_block_cycles(L) for L in neg]65    if not blocks:66        return 167    acc = blocks[0]68    for b in blocks[1:]:69        acc = _combine(acc, b)70    return sum(acc.values())7172def _partitions(n, m=None):73    if m is None:74        m = n75    if n == 0:76        yield ()77        return78    for k in range(min(n, m), 0, -1):79        for rest in _partitions(n - k, k):80            yield (k,) + rest8182def _bipartitions(dimension):83    for s in range(dimension + 1):84        for pos in (_partitions(s) if s > 0 else [()]):85            for neg in (_partitions(dimension - s) if (dimension - s) > 0 else [()]):86                yield pos, neg8788def _class_size(pos, neg, dimension):89    centralizer = 190    for part in (pos, neg):91        for length, mult in Counter(part).items():92            centralizer *= factorial(mult) * ((2 * length) ** mult)93    return (2 ** dimension * factorial(dimension)) // centralizer9495def total_designs(dimension):96    return 2 ** (2 ** dimension)9798def distinct_designs(dimension):99    total = Fraction(0)100    order = 2 ** dimension * factorial(dimension)101    checked = 0102    for pos, neg in _bipartitions(dimension):103        size = _class_size(pos, neg, dimension)104        checked += size105        total += size * (2 ** _class_cycles(pos, neg))106    if checked != order:107        raise ValueError("class sizes do not sum to the group order.")108    if total.denominator != 1:109        raise ValueError("Burnside average is not an integer.")110    return int(total // order)111112def sequence(max_dimension):113    return [distinct_designs(d) for d in range(1, max_dimension + 1)]