"""Exact tools for wallpaper groups, shared by the check scripts in this folder.

Python 3 standard library only. Every number is an int or a Fraction; nothing is floating point.

An isometry of the plane is stored in LATTICE COORDINATES as an affine map g = (M, t):
    v  |->  M v + t,
with M a 2x2 integer matrix (the lattice basis is fixed by each script) and t a pair of Fractions.
Composition reads like function composition: (g * h)(v) = g(h(v)), so h acts first.

A wallpaper group G is stored as one coset representative (P, s) per element P of its point group;
G is then { (P, s + w) : w in Z^2 }, i.e. the lattice is Z^2 in the chosen basis.
For n >= 1 the lattice nZ^2 is normal in G (every P maps it to itself), so G/nZ^2 is a finite group
of order |P| * n^2: an element is (P, t mod n). A finite quotient is how non-generation is proved:
if a set S maps to a proper subgroup of G/nZ^2, then S cannot generate G.
"""
from fractions import Fraction as F
from itertools import combinations

I2 = ((1, 0), (0, 1))


# ---------------------------------------------------------------- 2x2 integer linear algebra
def mm(A, B):
    return ((A[0][0] * B[0][0] + A[0][1] * B[1][0], A[0][0] * B[0][1] + A[0][1] * B[1][1]),
            (A[1][0] * B[0][0] + A[1][1] * B[1][0], A[1][0] * B[0][1] + A[1][1] * B[1][1]))


def mv(A, v):
    return (A[0][0] * v[0] + A[0][1] * v[1], A[1][0] * v[0] + A[1][1] * v[1])


def det(A):
    return A[0][0] * A[1][1] - A[0][1] * A[1][0]


def inv_unimodular(A):
    """Inverse of an integer matrix with determinant +-1 (integer again)."""
    d = det(A)
    assert d in (1, -1), f"not unimodular: {A}"
    return ((d * A[1][1], -d * A[0][1]), (-d * A[1][0], d * A[0][0]))


def vadd(u, v):
    return (u[0] + v[0], u[1] + v[1])


def vsub(u, v):
    return (u[0] - v[0], u[1] - v[1])


def vneg(u):
    return (-u[0], -u[1])


def fr(v):
    return (F(v[0]), F(v[1]))


def is_int_vec(v):
    return all(F(x).denominator == 1 for x in v)


# ---------------------------------------------------------------- affine maps
class Aff:
    __slots__ = ("M", "t")

    def __init__(self, M, t=(0, 0)):
        self.M = tuple(tuple(int(x) for x in row) for row in M)
        self.t = fr(t)

    def __mul__(self, o):  # self after o
        return Aff(mm(self.M, o.M), vadd(mv(self.M, o.t), self.t))

    def inv(self):
        Mi = inv_unimodular(self.M)
        return Aff(Mi, vneg(mv(Mi, self.t)))

    def key(self):
        return (self.M, self.t)

    def __eq__(self, o):
        return self.key() == o.key()

    def __hash__(self):
        return hash(self.key())

    def __call__(self, v):
        return vadd(mv(self.M, fr(v)), self.t)

    def __pow__(self, k):
        r = Aff(I2)
        base = self if k >= 0 else self.inv()
        for _ in range(abs(k)):
            r = r * base
        return r

    def mod(self, n):
        return (self.M, (self.t[0] % n, self.t[1] % n))

    def __repr__(self):
        def f(x):
            return str(x)
        return f"[M={self.M}, t=({f(self.t[0])},{f(self.t[1])})]"


def T(w):
    """Translation by w."""
    return Aff(I2, w)


def rot_about(R, c):
    """The isometry with linear part R that fixes the point c:  v -> R(v - c) + c."""
    c = fr(c)
    return Aff(R, vsub(c, mv(R, c)))


def fixed_point(g):
    """Unique fixed point of g = (M, t) when M - I is invertible (a rotation), as Fractions."""
    (a, b), (c, d) = g.M
    # solve (M - I) p = -t
    A = ((a - 1, b), (c, d - 1))
    D = det(A)
    assert D != 0
    rhs = vneg(g.t)
    return (F(A[1][1] * rhs[0] - A[0][1] * rhs[1], D), F(-A[1][0] * rhs[0] + A[0][0] * rhs[1], D))


# ---------------------------------------------------------------- words
def word_str(word, names=None):
    """('y',1),('x',-1)  ->  'y x^-1'."""
    if not word:
        return "1"
    out = []
    for n, e in word:
        out.append(n if e == 1 else f"{n}^-1")
    return " ".join(out)


def evaluate(word, gens):
    g = Aff(I2)
    for n, e in word:
        g = g * (gens[n] if e == 1 else gens[n].inv())
    return g


def find_words(gens, targets, max_nodes=400000):
    """Breadth-first search over words in gens and their inverses; returns {label: shortest word}
    for every target reached. Each word is re-checked by evaluate() by the caller."""
    letters = [(n, 1, g) for n, g in gens.items()] + [(n, -1, g.inv()) for n, g in gens.items()]
    start = Aff(I2)
    seen = {start.key(): ()}
    tkeys = {}
    for lab, tg in targets.items():
        tkeys.setdefault(tg.key(), []).append(lab)
    found = {}
    for lab in tkeys.get(start.key(), []):
        found[lab] = ()
    frontier = [start]
    while frontier and len(found) < len(targets) and len(seen) < max_nodes:
        nxt = []
        for g in frontier:
            w = seen[g.key()]
            for (n, e, a) in letters:
                h = g * a
                k = h.key()
                if k in seen:
                    continue
                seen[k] = w + ((n, e),)
                for lab in tkeys.get(k, []):
                    if lab not in found:
                        found[lab] = seen[k]
                nxt.append(h)
        frontier = nxt
    return found


# ---------------------------------------------------------------- wallpaper groups
class WGroup:
    def __init__(self, orb, iuc, reps, note=""):
        self.orb = orb
        self.iuc = iuc
        self.reps = list(reps)  # one Aff per point-group element
        self.note = note
        self.by_M = {r.M: r for r in self.reps}
        assert len(self.by_M) == len(self.reps), "two reps with the same linear part"

    def contains(self, g):
        r = self.by_M.get(g.M)
        return r is not None and is_int_vec(vsub(g.t, r.t))

    def check_is_group(self):
        """The reps' products and inverses land in G, the identity coset is the lattice, every P is
        in GL(2,Z): so G = union of (P, s + Z^2) is a group whose translation subgroup is Z^2."""
        assert I2 in self.by_M and is_int_vec(self.by_M[I2].t)
        for r in self.reps:
            assert det(r.M) in (1, -1)
            assert self.contains(r.inv()), (self.orb, r)
            for q in self.reps:
                assert self.contains(r * q), (self.orb, r, q)
        return True

    def standard_generators(self):
        """Generators of G by construction: the two lattice translations and one rep per coset."""
        out = {"t(1,0)": T((1, 0)), "t(0,1)": T((0, 1))}
        for i, r in enumerate(self.reps):
            if r.M != I2:
                out[f"rep{i}"] = r
        return out

    def order_mod(self, n):
        return len(self.reps) * n * n

    def elements_mod(self, n):
        out = []
        for r in self.reps:
            for i in range(n):
                for j in range(n):
                    out.append((r.M, ((r.t[0] + i) % n, (r.t[1] + j) % n)))
        assert len(set(out)) == self.order_mod(n)
        return out


def mul_mod(a, b, n):
    M = mm(a[0], b[0])
    t = vadd(mv(a[0], b[1]), a[1])
    return (M, (t[0] % n, t[1] % n))


def closure_mod(gen_keys, n):
    """The subgroup of G/nZ^2 generated by the given mod-n keys (finite, so right multiplication by
    the generators alone reaches everything)."""
    e = (I2, (F(0), F(0)))
    seen = {e}
    frontier = [e]
    while frontier:
        nxt = []
        for x in frontier:
            for g in gen_keys:
                y = mul_mod(x, g, n)
                if y not in seen:
                    seen.add(y)
                    nxt.append(y)
        frontier = nxt
    return seen


class FiniteQuotient:
    """G/nZ^2 with a multiplication table on indices, for exhaustive searches."""

    def __init__(self, group, n):
        self.n = n
        self.E = group.elements_mod(n)
        self.idx = {x: i for i, x in enumerate(self.E)}
        N = len(self.E)
        self.N = N
        self.tab = [[self.idx[mul_mod(self.E[i], self.E[j], n)] for j in range(N)] for i in range(N)]
        self.e = self.idx[(I2, (F(0), F(0)))]

    def generated(self, gens):
        seen = {self.e}
        frontier = [self.e]
        tab = self.tab
        while frontier:
            nxt = []
            for x in frontier:
                row = tab[x]
                for g in gens:
                    y = row[g]
                    if y not in seen:
                        seen.add(y)
                        nxt.append(y)
            frontier = nxt
        return seen

    def some_k_generate(self, k):
        """True if some k elements generate the whole finite group (exhaustive)."""
        for combo in combinations(range(self.N), k):
            if len(self.generated(combo)) == self.N:
                return True
        return False

    def is_abelian(self):
        t = self.tab
        return all(t[i][j] == t[j][i] for i in range(self.N) for j in range(self.N))

    def exponent_divides(self, m):
        for i in range(self.N):
            x = self.e
            for _ in range(m):
                x = self.tab[x][i]
            if x != self.e:
                return False
        return True


# ---------------------------------------------------------------- Conway's Magic Theorem
def magic_cost(symbol):
    """Cost of an orbifold symbol: o = 2, * = 1, x = 1, a gyration n = (n-1)/n, a corner n = (n-1)/(2n).
    Digits before any * or x are gyrations; digits after a * are corners."""
    cost = F(0)
    seen_star = False
    for ch in symbol:
        if ch == "o":
            cost += 2
        elif ch == "*":
            cost += 1
            seen_star = True
        elif ch == "x":
            cost += 1
        elif ch.isdigit():
            n = int(ch)
            cost += F(n - 1, 2 * n) if seen_star else F(n - 1, n)
        else:
            raise ValueError(ch)
    return cost


def fmt_fr(x):
    x = F(x)
    return str(x.numerator) if x.denominator == 1 else f"{x.numerator}/{x.denominator}"