"""Rigid triangles and commutative relations (Sections 6 and 8).
A homogeneous coherent configuration is given by a matrix R, with R[i][k] being the relation of the pair (i,k).
  Lemma 6.2: vanishing of the products u_ij u_kl u_ij', in 2 situations (coherent versions of 2 lemmas of Schmidt).
  Lemma 6.3: a relation is commutative when its valency is 1, when its pairs are part of rigid triples, and so on.
If all the relations are commutative, there is no quantum symmetry (Theorem 6.4).  These lemmas are coherent algebra
versions of results of Schmidt, and the closure computed here uses Lemma 6.3 (1) and (4).  The script checks that this is the
case under the bounds of Theorems 8.1 and 8.2, for the cyclotomic schemes on Z_N.   Output: rigid.log"""
import sys
from schur import *

def cayley_relations(N, sub, col):
    return [[col[sub[y][x]] for y in range(N)] for x in range(N)]

def commutative_relations(R, start=()):
    """Closure of the set 'start' of relations under Lemmas 6.2 and 6.3.  The configuration is assumed to be
    vertex-transitive, so that each relation consists of the translates of its pairs (0,l).
    For a pair (j,l) = (0,l) the products u_ij u_kl u_ij', with j' != j, must all vanish.  This is the case when
    there is a vertex n with R[n][j'] != R[n][j], such that the relation of (l,n) is commutative (Lemma 6.2 (1)),
    or such that l is the only vertex having its relations to j, j', n (Lemma 6.2 (2))."""
    N = len(R)
    nb = [{} for _ in range(N)]
    for i in range(N):
        for v in range(N): nb[i].setdefault(R[v][i], []).append(v)      # nb[i][r] = {v : R[v][i] = r}
    tr = {R[0][k]: R[k][0] for k in range(N)}
    cls = {}
    for l in range(1, N): cls.setdefault(R[0][l], []).append(l)
    comm = set(start) | {R[0][0]}
    for r, ls in cls.items():
        if len(ls) == 1: comm.add(r)
    comm |= {tr[r] for r in comm}
    def killed(l, p):                    # j = 0, j' = p: is u_i0 u_kl u_ip = 0 ?
        cand = nb[0][R[l][0]]
        for n in range(N):
            if R[n][p] == R[n][0]: continue
            if R[l][n] in comm: return True
            a, b = R[l][p], R[l][n]
            if sum(1 for t in cand if R[t][p] == a and R[t][n] == b) == 1: return True
        return False
    changed = True
    while changed:
        changed = False
        for r, ls in cls.items():
            if r in comm: continue
            if all(killed(l, p) for l in ls for p in range(1, N) if R[l][p] == R[l][0]):
                comm |= {r, tr[r]}; changed = True
    return comm

def all_commutative(R, start=()):
    return len(commutative_relations(R, start)) == len({r for row in R for r in row})

def orbit_coloring(N, E):
    col = [-1] * N; col[0] = 0; c = 0
    for x in range(1, N):
        if col[x] < 0:
            c += 1
            for a in E: col[a * x % N] = c
    return col

if __name__ == "__main__":
    out = open("rigid.log", "w")
    def say(*a):
        print(*a); print(*a, file=out); out.flush()
    say("DIHEDRAL TYPE on Z_N (E = {+-1}): all the relations commutative?")
    for N in range(3, 21):
        _, neg, sub, cl = setup((N,))
        say("  N=%2d: %s" % (N, all_commutative(cayley_relations(N, sub, orbit_coloring(N, [1, N - 1])))))
    say("PRIME CASE, p < 200: number of pairs (p,k) where all the relations are commutative,")
    say("  compared with the bound p > (k-1)(k-2)+2 of Theorem 8.2")
    tot = cov = bound = 0; viol = []
    for p in primes(200):
        if p < 5: continue
        _, neg, sub, cl = setup((p,))
        for k in range(2, p - 1, 2):
            if (p - 1) % k: continue
            r = all_commutative(cayley_relations(p, sub, orbit_coloring(p, cyclic_subgroup(p, k))))
            b = p > (k - 1) * (k - 2) + 2
            tot += 1; cov += r; bound += b
            if b and not r: viol.append((p, k))
    say("  total %d, commutative %d, bound applies %d, violations of the bound: %s" % (tot, cov, bound, viol))
    say("GENERAL CASE, N < 130, all the groups -1 in E generated by at most 3 elements, k = |E|, m = max gcd(a-1,N):")
    say("  all the relations must be commutative when N > mk(k-1)-2(k-2) (Theorem 8.2), and when N > 2(k-1)(km-1)+2 (Theorem 8.1)")
    groups = ok = b1 = b2 = 0; viol = []
    for N in range(5, 130):
        _, neg, sub, cl = setup((N,))
        units = [a for a in range(1, N) if math.gcd(a, N) == 1]
        done = set()
        for a in units:
            for b in units:
                if b < a: continue
                E = {1, N - 1}
                while True:
                    F = E | {x * g % N for x in E for g in (a, b)}
                    if F == E: break
                    E = F
                E = tuple(sorted(E))
                if E in done or len(E) > 12: continue
                done.add(E); groups += 1
                k = len(E); m = max(math.gcd(c - 1, N) for c in E if c != 1)
                r = all_commutative(cayley_relations(N, sub, orbit_coloring(N, E)))
                t2 = N > m * k * (k - 1) - 2 * (k - 2); t1 = N > 2 * (k - 1) * (k * m - 1) + 2
                ok += r; b2 += t2; b1 += t1
                if t2 and not r: viol.append((N, E))
    say("  groups %d, commutative %d, Theorem 8.2 applies %d, Theorem 8.1 applies %d, violations: %s" % (groups, ok, b2, b1, viol))
    say("PRIME POWERS, E of order k | p-1:")
    for N, p in [(9, 3), (25, 5), (27, 3), (49, 7), (121, 11), (125, 5), (169, 13)]:
        _, neg, sub, cl = setup((N,))
        for k in range(2, p, 2):
            if (p - 1) % k: continue
            say("  N=%3d k=%2d: %s" % (N, k, all_commutative(cayley_relations(N, sub, orbit_coloring(N, cyclic_subgroup(N, k))))))
    out.close()