"""Basic tools: finite abelian groups, Cayley graphs, coherent closure (Schur ring generated by S),
automorphism counting, subgroups of Z_N^*.  Used by all the other scripts."""
import itertools, math

def phi(n):
    r, d, m = n, 2, n
    while d * d <= m:
        if m % d == 0:
            r = r // d * (d - 1)
            while m % d == 0: m //= d
        d += 1
    if m > 1: r = r // m * (m - 1)
    return r

def primes(n):
    return [q for q in range(2, n) if all(q % d for d in range(2, int(q ** 0.5) + 1))]

def cyclic_subgroup(N, k):
    """The subgroup of order k of Z_N^*, when Z_N^* is cyclic."""
    ph = phi(N)
    for g in range(2, N):
        if math.gcd(g, N) > 1: continue
        o, x = 1, g
        while x != 1: x = x * g % N; o += 1
        if o == ph:
            h = pow(g, ph // k, N); return sorted({pow(h, i, N) for i in range(k)})
    raise ValueError("Z_N^* is not cyclic")

def setup(dims):
    """A = Z_{d1} x ... x Z_{dr}.  Returns N, neg (inversion), sub[a][x] = a - x, and the classes {x,-x}, x != 0."""
    els = list(itertools.product(*[range(n) for n in dims]))
    idx = {x: i for i, x in enumerate(els)}
    N = len(els)
    neg = [idx[tuple((-a) % n for a, n in zip(x, dims))] for x in els]
    sub = [[idx[tuple((a - b) % n for a, b, n in zip(x, y, dims))] for y in els] for x in els]
    cl, seen = [], set()
    for i in range(1, N):
        if i in seen: continue
        cl.append((i, neg[i])); seen |= {i, neg[i]}
    return N, neg, sub, cl

def closure(N, sub, col):
    """Weisfeiler-Leman on a Cayley graph: refines the partition col of A until it is a Schur ring.
    Returns the final partition; its number of classes is the rank of the coherent algebra."""
    while True:
        k = max(col) + 1
        sig = []
        for a in range(N):
            cnt = [0] * (k * k)
            sa = sub[a]
            for x in range(N):
                cnt[col[x] * k + col[sa[x]]] += 1
            sig.append((col[a], tuple(cnt)))
        names = {s: i for i, s in enumerate(sorted(set(sig)))}
        new = [names[s] for s in sig]
        if max(new) + 1 == k: return new
        col = new

def coloring(N, S):
    col = [2] * N; col[0] = 0
    for s in S: col[s] = 1
    return col

def rank(N, sub, S):
    return max(closure(N, sub, coloring(N, S))) + 1

def aut_order(N, adj):
    """Order of the automorphism group of a vertex-transitive graph (N times the stabilizer of 0)."""
    img = [-1] * N; used = [False] * N
    def rec(i):
        if i == N: return 1
        tot = 0
        for w in range(N):
            if used[w]: continue
            if all(adj[i][j] == adj[w][img[j]] for j in range(i)):
                img[i] = w; used[w] = True; tot += rec(i + 1); used[w] = False
        img[i] = -1
        return tot
    img[0] = 0; used[0] = True
    return N * rec(1)

def sets_from_mask(mask, cl):
    S = set()
    for c in range(len(cl)):
        if mask >> c & 1: S |= set(cl[c])
    return S