#!/usr/bin/env python3
"""Finite permutation groups, with the three places they earn their keep in a plant.

A permutation of n things is a bijection {0..n-1} -> {0..n-1}, stored as a tuple p with p[i] = image of i.
Composition (p * q)(i) = p[q[i]]: apply q first. A permutation group is a set of permutations closed under
composition and inverse. Everything below is built from that and checked by brute force.

  1. Orders and cycles: a perfect riffle shuffle of 52 cards is a permutation; its order (lcm of its cycle
     lengths) is how many shuffles return the deck. Out-shuffle: 8. In-shuffle: 52.
  2. Burnside's lemma: the number of distinct colorings up to a symmetry group G is
         (1/|G|) * sum over g in G of  k^(number of cycles of g).
     Box panels under rotation (C4) and rotation+flip (D4); die positions on a sheet under the sheet's
     symmetries (the Klein four-group V).
  3. Verhoeff's check digit: the dihedral group D5 (order 10) used as the "addition" of a checksum.
     It catches every single-digit error and every adjacent transposition, which no mod-10 scheme can.

Run:  python3 perms.py          # prints the numbers the post uses and writes perms.json
      python3 perms.py --test   # brute-force checks of every claim
"""
from __future__ import annotations

import itertools
import json
import math
import random
import sys
from pathlib import Path

HERE = Path(__file__).resolve().parent


# ----------------------------------------------------------------------------- permutations

def compose(p, q):
    """(p * q)(i) = p[q[i]]: do q, then p."""
    return tuple(p[q[i]] for i in range(len(p)))


def inverse(p):
    inv = [0] * len(p)
    for i, j in enumerate(p):
        inv[j] = i
    return tuple(inv)


def identity(n):
    return tuple(range(n))


def cycles(p):
    """Cycle decomposition, fixed points included as 1-cycles."""
    seen, out = set(), []
    for i in range(len(p)):
        if i in seen:
            continue
        c, j = [], i
        while j not in seen:
            seen.add(j)
            c.append(j)
            j = p[j]
        out.append(tuple(c))
    return out


def order(p):
    o = 1
    for c in cycles(p):
        o = o * len(c) // math.gcd(o, len(c))
    return o


def generate(gens):
    """Closure of a set of generators under composition: the group they generate."""
    n = len(gens[0])
    group = {identity(n)}
    frontier = [identity(n)]
    while frontier:
        nxt = []
        for g in frontier:
            for h in gens:
                for x in (compose(g, h), compose(h, g)):
                    if x not in group:
                        group.add(x)
                        nxt.append(x)
        frontier = nxt
    return sorted(group)


def cycle_notation(p):
    parts = ["(" + " ".join(str(i + 1) for i in c) + ")" for c in cycles(p) if len(c) > 1]
    return "".join(parts) or "e"


# ----------------------------------------------------------------------------- 1. shuffles

def out_shuffle(n):
    """Cut a deck of n (even) in half and interleave, top card staying on top. Position i -> ?"""
    half = n // 2
    p = [0] * n
    for i in range(half):
        p[i] = 2 * i          # top half goes to even positions
        p[half + i] = 2 * i + 1
    return tuple(p)


def in_shuffle(n):
    half = n // 2
    p = [0] * n
    for i in range(half):
        p[i] = 2 * i + 1      # top half goes to odd positions: the top card moves down one
        p[half + i] = 2 * i
    return tuple(p)


# ----------------------------------------------------------------------------- 2. symmetry groups and Burnside

def cyclic_group(n):
    """C_n acting on n positions in a ring: generated by the rotation i -> i+1."""
    return generate([tuple((i + 1) % n for i in range(n))])


def dihedral_group(n):
    """D_n on n positions in a ring: rotation and a reflection."""
    rot = tuple((i + 1) % n for i in range(n))
    ref = tuple((-i) % n for i in range(n))
    return generate([rot, ref])


def rectangle_group(rows, cols):
    """The symmetries of an r x c sheet of die positions (r != c): identity, 180-degree turn, two flips. Klein four."""
    idx = lambda r, c: r * cols + c
    n = rows * cols
    turn = tuple(idx(rows - 1 - r, cols - 1 - c) for r in range(rows) for c in range(cols))
    flip_h = tuple(idx(rows - 1 - r, c) for r in range(rows) for c in range(cols))
    return generate([turn, flip_h]), {"turn 180": turn, "flip top-bottom": flip_h}


def burnside(group, k):
    """Number of orbits of k-colorings of the positions under the group: average number of fixed colorings."""
    total = sum(k ** len(cycles(g)) for g in group)
    assert total % len(group) == 0
    return total // len(group)


def orbits(group, k):
    """Brute force: partition all k^n colorings into orbits. For checking Burnside and for the figure."""
    n = len(group[0])
    seen, out = set(), []
    for col in itertools.product(range(k), repeat=n):
        if col in seen:
            continue
        orb = {tuple(col[g[i]] for i in range(n)) for g in group}
        seen |= orb
        out.append(sorted(orb))
    return out


# ----------------------------------------------------------------------------- 3. Verhoeff

# D5, the symmetries of a regular pentagon, as a 10 x 10 multiplication table. Elements 0-4 are the rotations
# (0 the identity), 5-9 the reflections. This is the table from Verhoeff's 1969 paper.
D5 = [
    [0, 1, 2, 3, 4, 5, 6, 7, 8, 9],
    [1, 2, 3, 4, 0, 6, 7, 8, 9, 5],
    [2, 3, 4, 0, 1, 7, 8, 9, 5, 6],
    [3, 4, 0, 1, 2, 8, 9, 5, 6, 7],
    [4, 0, 1, 2, 3, 9, 5, 6, 7, 8],
    [5, 9, 8, 7, 6, 0, 4, 3, 2, 1],
    [6, 5, 9, 8, 7, 1, 0, 4, 3, 2],
    [7, 6, 5, 9, 8, 2, 1, 0, 4, 3],
    [8, 7, 6, 5, 9, 3, 2, 1, 0, 4],
    [9, 8, 7, 6, 5, 4, 3, 2, 1, 0],
]
# The fixed permutation p applied to digits by position, and its powers. Chosen so that for every pair x != y,
# x * p(y) != y * p(x) in D5: that is exactly what makes adjacent transpositions detectable.
P = [[0, 1, 2, 3, 4, 5, 6, 7, 8, 9],
     [1, 5, 7, 6, 2, 8, 3, 0, 9, 4],
     [5, 8, 0, 3, 7, 9, 6, 1, 4, 2],
     [8, 9, 1, 6, 0, 4, 3, 5, 2, 7],
     [9, 4, 5, 3, 1, 2, 6, 8, 7, 0],
     [4, 2, 8, 6, 5, 7, 3, 9, 0, 1],
     [2, 7, 9, 3, 8, 0, 6, 4, 1, 5],
     [7, 0, 4, 6, 9, 1, 3, 2, 5, 8]]
INV = [0, 4, 3, 2, 1, 5, 6, 7, 8, 9]   # inverses in D5


def verhoeff_check(digits: str) -> int:
    """The check digit for a string of digits."""
    c = 0
    for i, d in enumerate(reversed(digits)):
        c = D5[c][P[(i + 1) % 8][int(d)]]
    return INV[c]


def verhoeff_valid(number_with_check: str) -> bool:
    c = 0
    for i, d in enumerate(reversed(number_with_check)):
        c = D5[c][P[i % 8][int(d)]]
    return c == 0


def luhn_check(digits: str) -> int:
    s = 0
    for i, d in enumerate(reversed(digits)):
        v = int(d)
        if i % 2 == 0:
            v *= 2
            if v > 9:
                v -= 9
        s += v
    return (10 - s % 10) % 10


def luhn_valid(number_with_check: str) -> bool:
    s = 0
    for i, d in enumerate(reversed(number_with_check)):
        v = int(d)
        if i % 2 == 1:
            v *= 2
            if v > 9:
                v -= 9
        s += v
    return s % 10 == 0


def sum_check(digits: str) -> int:
    return (10 - sum(int(d) for d in digits) % 10) % 10


def sum_valid(number_with_check: str) -> bool:
    return sum(int(d) for d in number_with_check) % 10 == 0


SCHEMES = {"Verhoeff (D5)": (verhoeff_check, verhoeff_valid), "Luhn (mod 10, doubling)": (luhn_check, luhn_valid), "plain digit sum mod 10": (sum_check, sum_valid)}


def error_catalog(number: str):
    """Every single-digit substitution and every adjacent transposition of a number (with its check digit)."""
    subs = [number[:i] + str(d) + number[i + 1:] for i in range(len(number)) for d in range(10) if str(d) != number[i]]
    trans = [number[:i] + number[i + 1] + number[i] + number[i + 2:] for i in range(len(number) - 1) if number[i] != number[i + 1]]
    return subs, trans


def detection_rates(sample=2000, length=7, seed=4):
    rng = random.Random(seed)
    out = {}
    for name, (check, valid) in SCHEMES.items():
        s_caught = s_total = t_caught = t_total = 0
        for _ in range(sample):
            base = "".join(rng.choice("0123456789") for _ in range(length))
            full = base + str(check(base))
            subs, trans = error_catalog(full)
            s_total += len(subs); t_total += len(trans)
            s_caught += sum(1 for e in subs if not valid(e))
            t_caught += sum(1 for e in trans if not valid(e))
        out[name] = {"single": s_caught / s_total, "adjacent": t_caught / t_total, "n_single": s_total, "n_adjacent": t_total}
    return out


# ----------------------------------------------------------------------------- run

def main(argv):
    C4, D4 = cyclic_group(4), dihedral_group(4)
    V, Vgens = rectangle_group(2, 3)
    if "--test" in argv:
        # group axioms by brute force
        for G in (C4, D4, V, dihedral_group(5)):
            S = set(G)
            assert all(compose(a, b) in S for a in G for b in G) and all(inverse(a) in S for a in G)
        assert len(C4) == 4 and len(D4) == 8 and len(V) == 4 and len(dihedral_group(5)) == 10
        # Burnside equals brute-force orbit counts
        for G in (C4, D4, V):
            for k in (2, 3):
                assert burnside(G, k) == len(orbits(G, k)), (G, k)
        assert burnside(C4, 2) == 6 and burnside(D4, 2) == 6 and burnside(C4, 3) == 24 and burnside(D4, 3) == 21 and burnside(V, 2) == 24
        # shuffles
        assert order(out_shuffle(52)) == 8 and order(in_shuffle(52)) == 52
        d = list(range(52)); p = out_shuffle(52)
        for _ in range(8):
            d = [d[inverse(p)[i]] for i in range(52)]
        assert d == list(range(52))
        # D5 is a group with the given table; P's defining property; Verhoeff catches everything it should
        assert all(sorted(row) == list(range(10)) for row in D5) and all(D5[a][INV[a]] == 0 for a in range(10))
        for a in range(10):
            for b in range(10):
                for c in range(10):
                    assert D5[D5[a][b]][c] == D5[a][D5[b][c]]
        assert all(D5[x][P[1][y]] != D5[y][P[1][x]] for x in range(10) for y in range(10) if x != y), "P's defining property"
        for i in range(1, 8):
            assert P[i] == tuple(P[1][x] for x in P[i - 1]) or list(P[i]) == [P[1][x] for x in P[i - 1]]
        rng = random.Random(1)
        for _ in range(300):
            base = "".join(rng.choice("0123456789") for _ in range(rng.randint(3, 10)))
            full = base + str(verhoeff_check(base))
            assert verhoeff_valid(full)
            subs, trans = error_catalog(full)
            assert all(not verhoeff_valid(e) for e in subs), "a single-digit error got through"
            assert all(not verhoeff_valid(e) for e in trans), "an adjacent transposition got through"
            lfull = base + str(luhn_check(base)); assert luhn_valid(lfull)
            sfull = base + str(sum_check(base)); assert sum_valid(sfull)
        assert luhn_valid("90" + str(luhn_check("09"))) or luhn_valid("09" + str(luhn_check("09")))  # Luhn's famous blind spot: 09 <-> 90
        rates = detection_rates(sample=300)
        assert rates["Verhoeff (D5)"]["single"] == 1.0 and rates["Verhoeff (D5)"]["adjacent"] == 1.0
        assert rates["plain digit sum mod 10"]["adjacent"] == 0.0
        print("groups close; Burnside == brute-force orbits; out-shuffle order 8, in-shuffle 52; D5 associative; Verhoeff catches 100% of single errors and adjacent transpositions; digit-sum catches 0% of transpositions")
        return 0

    rates = detection_rates()
    out = {
        "C4": [list(g) for g in C4], "D4": [list(g) for g in D4], "V23": [list(g) for g in V],
        "V23_generators": {k: list(v) for k, v in Vgens.items()},
        "burnside": {"C4": {k: burnside(C4, k) for k in (2, 3, 4)}, "D4": {k: burnside(D4, k) for k in (2, 3, 4)}, "V23": {k: burnside(V, k) for k in (2, 3)}},
        "shuffles": {str(n): {"out": {"order": order(out_shuffle(n)), "cycles": [len(c) for c in cycles(out_shuffle(n))]}, "in": {"order": order(in_shuffle(n)), "cycles": [len(c) for c in cycles(in_shuffle(n))]}} for n in (8, 10, 24, 52)},
        "D5": D5, "P": P, "INV": INV,
        "rates": rates,
    }
    (HERE / "perms.json").write_text(json.dumps(out, separators=(",", ":")))
    print("Burnside counts:")
    print(f"  box panels, rotation only (C4):      2 colors {burnside(C4,2)}, 3 colors {burnside(C4,3)}, 4 colors {burnside(C4,4)}")
    print(f"  box panels, rotation + flip (D4):    2 colors {burnside(D4,2)}, 3 colors {burnside(D4,3)}, 4 colors {burnside(D4,4)}")
    print(f"  2x3 sheet of die positions (V):      blank/waste {burnside(V,2)} of 64 layouts are distinct")
    for n in (8, 24, 52):
        print(f"  deck of {n}: out-shuffle order {order(out_shuffle(n))} (cycles {sorted(len(c) for c in cycles(out_shuffle(n)))}), in-shuffle order {order(in_shuffle(n))}")
    print("detection rates over 2,000 random 7-digit numbers with check digit:")
    for name, r in rates.items():
        print(f"  {name:<26} single-digit errors {r['single']:6.1%}   adjacent transpositions {r['adjacent']:6.1%}")
    print(f"wrote {HERE / 'perms.json'}")
    return 0


if __name__ == "__main__":
    raise SystemExit(main(sys.argv[1:]))