"""Indistinguishing numbers and pseudocyclic schemes (Theorem 8.3).
For a homogeneous coherent configuration R, with maximal valency k, let c be the maximal number of vertices v
satisfying R[i][v] = R[i'][v], over the pairs of distinct vertices i, i'.  The proof of Theorem 8.3 shows that
all the relations are commutative, by Lemma 6.3 (2), when N > 2(k-1)(c-1)+2.  The script
  (1) checks that c = k-1 for the cyclotomic schemes over Z_p (Frobenius schemes, which are pseudocyclic),
      including those with -1 not in E, which are not symmetric, and tests the bound N > 2(k-1)(k-2)+2,
  (2) tests the bound N > 2(k-1)(c-1)+2 for the cyclotomic schemes over Z_N, and compares c with km,
  (3) does the same for the coset schemes of transitive.py,
  (4) checks the proofs of Theorems 8.1 and 8.2 directly: under the bounds there, any pair of distinct vertices (i,j)
      has a vertex v such that i is determined by its relations with j,v, and j by its relations with i,v.
Output: pseudocyclic.log"""
import math
from schur import *
from rigid import cayley_relations, commutative_relations, all_commutative, orbit_coloring
from transitive import coset_action, orbitals, perm, mul

def indist(R, transitive=True):
    """The number c.  For a vertex-transitive configuration it is enough to take i = 0."""
    N = len(R)
    return max(sum(1 for v in range(N) if R[i][v] == R[j][v]) for i in ([0] if transitive else range(N)) for j in range(N) if j != i)

def direct(R):
    """Does Lemma 6.3 (2) apply to all the pairs (0,l)?"""
    N = len(R)
    def det(i, k, m): return sum(1 for t in range(N) if R[t][k] == R[i][k] and R[t][m] == R[i][m]) == 1
    return all(any(det(0, l, v) and det(l, 0, v) for v in range(1, N) if v != l) for l in range(1, N))

def max_valency(R):
    cnt = {}
    for v in range(1, len(R)): cnt[R[0][v]] = cnt.get(R[0][v], 0) + 1
    return max(cnt.values())

if __name__ == "__main__":
    out = open("pseudocyclic.log", "w")
    def say(*a):
        print(*a); print(*a, file=out); out.flush()
    say("(1) CYCLOTOMIC SCHEMES OVER Z_p, p < 150, all the groups E of order 2 <= k < p-1:")
    tot = sym = eq = bound = comm = 0; viol = []
    for p in primes(150):
        if p < 5: continue
        _, neg, sub, cl = setup((p,))
        for k in range(2, p - 1):
            if (p - 1) % k: continue
            E = cyclic_subgroup(p, k)
            R = cayley_relations(p, sub, orbit_coloring(p, E))
            c = indist(R); r = all_commutative(R); b = p > 2 * (k - 1) * (k - 2) + 2
            tot += 1; sym += (p - 1) in E; eq += c == k - 1; bound += b; comm += r
            if b and not r: viol.append((p, k))
    say("  schemes %d (symmetric %d), with c = k-1: %d, bound N > 2(k-1)(k-2)+2 applies %d, all commutative %d, violations: %s" % (tot, sym, eq, bound, comm, viol))
    say("(2) CYCLOTOMIC SCHEMES OVER Z_N, N < 100, groups -1 in E generated by at most 3 elements:")
    groups = le = bound = comm = b82 = d82 = 0; viol = []
    for N in range(5, 100):
        _, 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(x - 1, N) for x in E if x != 1)
                R = cayley_relations(N, sub, orbit_coloring(N, E))
                c = indist(R); kk = max_valency(R); r = all_commutative(R); t = N > 2 * (kk - 1) * (c - 1) + 2
                le += c <= k * m; bound += t; comm += r
                if N > m * k * (k - 1) - 2 * (k - 2): b82 += 1; d82 += direct(R)
                if t and not r: viol.append((N, E))
    say("  groups %d, with c <= km: %d, bound N > 2(k-1)(c-1)+2 applies %d, all commutative %d, violations: %s" % (groups, le, bound, comm, viol))
    say("  (4): bound of Theorem 8.2 applies %d, Lemma 6.3 (2) applies to all the pairs %d" % (b82, d82))
    say("(3) COSET SCHEMES, with (4):")
    t7 = tuple([(x + 1) % 7 for x in range(7)] + [7])
    inv = {0: 7, 7: 0}
    for x in range(1, 7): inv[x] = (-pow(x, 5, 7)) % 7
    s7 = tuple(inv[x] for x in range(8))
    a5 = [perm(5, (0, 1, 2, 3, 4)), perm(5, (0, 1, 2))]
    s5 = [perm(5, (0, 1, 2, 3, 4)), perm(5, (0, 1))]
    for name, gens, hgens in [("PSL(2,7) / Z_2", [t7, s7], [s7]), ("PSL(2,7) / Z_3", [t7, s7], [mul(s7, t7)]),
                              ("A_5 / Z_2", a5, [perm(5, (0, 1), (2, 3))]), ("A_5 / Z_3", a5, [perm(5, (0, 1, 2))]),
                              ("S_5 / <(12)(34)>", s5, [perm(5, (0, 1), (2, 3))]), ("S_5 / <(123)>", s5, [perm(5, (0, 1, 2))])]:
        G, H, N, act = coset_action(gens, hgens)
        k = len(H); e = tuple(range(len(gens[0])))
        m = max(sum(1 for i in range(N) if act[g][i] == i) for g in G if g != e)
        R, rank = orbitals(N, act, gens)
        c = indist(R); kk = max_valency(R); b = 2 * (kk - 1) * (c - 1) + 2
        r = len(commutative_relations(R)) == rank
        assert r or N <= b
        d = direct(R)
        assert d or N <= 2 * (k - 1) * (k * m - 1) + 2
        say("  %-18s N=%3d k=%d m=%d c=%2d (km=%2d)  bound=%3d applies=%-5s all commutative=%-5s Lemma 6.3 (2) for all the pairs=%s" % (name, N, k, m, c, k * m, b, N > b, r, d))
    out.close()