#!/usr/bin/env python3
"""Ancillary verification script for the note
"Hadamard difference sets in groups with a cyclic odd quotient"
(David Niedbala Giraudin).

Exact integer arithmetic, except in Part B, where a floating-point filter on the
character values (tolerance 1e-6, while the values have absolute value 3) only
discards candidates; every survivor is then checked exactly.  Requires numpy.
All parameters are fixed.  Expected running time: about one minute.

Parts
  0. self-test of the cyclotomic arithmetic
  A. the algebraic identities of the Klein and cyclic 2-parts on random sign arrays
  B. exhaustive check of the alternating-product proposition for all x in
     Z[C_9] and all x in Z[i][C_9] with x x* = 9 (integrality, root of unity,
     3-power order, residue 1 above 3)
  C. two explicit (36,15,6) difference sets: where the argument breaks when the
     hypothesis "cyclic Sylow subgroup of full order" is violated
  D. an explicit element of Z[i][C_3] showing why self-conjugacy is assumed
  E. counts of abelian groups of order 4u**2 (u odd, 3 <= u <= 105) covered by
     the main theorem
  F. the values u <= 200000 for which C2 x C_{2u**2} is not already excluded by
     two classical criteria (exponent bound under self-conjugacy; field descent)
"""
import itertools
import random
import sys
import time

T0 = time.time()


def say(*a):
    print(*a, flush=True)


# ---------------------------------------------------------------- arithmetic
def poly_divmod(a, b):
    """Integer polynomial division by a monic b (coefficient lists, low first)."""
    a = list(a)
    db = len(b) - 1
    assert b[-1] == 1
    q = [0] * max(len(a) - db, 1)
    for k in range(len(a) - 1, db - 1, -1):
        c = a[k]
        if c:
            q[k - db] = c
            for j in range(db + 1):
                a[k - db + j] -= c * b[j]
    r = a[:db] if db > 0 else []
    while len(r) < db:
        r.append(0)
    return q, r


_PHI = {}


def cyclotomic(n):
    if n in _PHI:
        return _PHI[n]
    num = [-1] + [0] * (n - 1) + [1]  # X**n - 1
    for d in range(1, n):
        if n % d == 0:
            q, r = poly_divmod(num, cyclotomic(d))
            assert all(c == 0 for c in r)
            num = q
    while num and num[-1] == 0:
        num.pop()
    _PHI[n] = num
    return num


def gcd(a, b):
    while b:
        a, b = b, a % b
    return abs(a)


class Cyc:
    """Elements of Z[zeta_N] in the power basis 1, zeta, ..., zeta**(phi(N)-1)."""

    def __init__(self, N):
        self.N = N
        self.phi = cyclotomic(N)
        self.deg = len(self.phi) - 1
        self.units = [k for k in range(1, N + 1) if gcd(k, N) == 1]

    def from_exps(self, terms):
        """terms: list of (integer coefficient, exponent of zeta)."""
        v = [0] * self.N
        for c, e in terms:
            v[e % self.N] += c
        return self.reduce(v)

    def reduce(self, v):
        if len(v) <= self.deg:
            return list(v) + [0] * (self.deg - len(v))
        return poly_divmod(v, self.phi)[1]

    def mul(self, a, b):
        v = [0] * (len(a) + len(b))
        for i, x in enumerate(a):
            if x:
                for j, y in enumerate(b):
                    if y:
                        v[i + j] += x * y
        return self.reduce(v)

    def add(self, a, b):
        return [x + y for x, y in zip(a, b)]

    def sub(self, a, b):
        return [x - y for x, y in zip(a, b)]

    def one(self):
        return self.from_exps([(1, 0)])

    def zeta(self, e):
        return self.from_exps([(1, e)])

    def sigma(self, a, k):
        return self.from_exps([(c, i * k) for i, c in enumerate(a) if c])

    def conj(self, a):
        return self.sigma(a, self.N - 1)

    def norm(self, a):
        r = self.one()
        for k in self.units:
            r = self.mul(r, self.sigma(a, k))
        assert all(c == 0 for c in r[1:])
        return r[0]

    def divide(self, a, b):
        """Return a/b if it lies in Z[zeta_N], else None (exact)."""
        nb = self.norm(b)
        assert nb != 0
        cof = self.one()
        for k in self.units:
            if k != 1:
                cof = self.mul(cof, self.sigma(b, k))
        num = self.mul(a, cof)
        if all(c % nb == 0 for c in num):
            return [c // nb for c in num]
        return None

    def is_zero(self, a):
        return all(c == 0 for c in a)

    def root_order(self, a):
        """Order of a if a is a root of unity of Q(zeta_N), else None."""
        cands = []
        for j in range(self.N):
            cands.append((self.zeta(j), self.N // gcd(j, self.N)))
            if self.N % 2 == 1:
                negz = [-c for c in self.zeta(j)]
                o = self.N // gcd(j, self.N)
                cands.append((negz, 2 * o))
        for z, o in cands:
            if z == a:
                return o
        return None


def rad(n):
    r, p = 1, 2
    while p * p <= n:
        if n % p == 0:
            r *= p
            while n % p == 0:
                n //= p
        p += 1
    if n > 1:
        r *= n
    return r


def factor(n):
    f, p = {}, 2
    while p * p <= n:
        while n % p == 0:
            f[p] = f.get(p, 0) + 1
            n //= p
        p += 1
    if n > 1:
        f[n] = f.get(n, 0) + 1
    return f


# group-ring elements of Z[i][C_M]: dict exponent -> (re, im)
def gr_mul(x, y, M):
    r = {}
    for a, (p, q) in x.items():
        for b, (s, t) in y.items():
            k = (a + b) % M
            u, v = r.get(k, (0, 0))
            r[k] = (u + p * s - q * t, v + p * t + q * s)
    return {k: v for k, v in r.items() if v != (0, 0)}


def gr_star(x, M):
    return {(-a) % M: (p, -q) for a, (p, q) in x.items()}


def char_value(x, M, primes, S, K):
    """chi_S(x) in Z[zeta_N]; chi_S(X) = zeta_M ** (sum_{p in S} M/p)."""
    N = K.N
    a = sum(M // p for p in S)
    terms = []
    for z, (re, im) in x.items():
        e = (z * a % M) * (N // M) if N % M == 0 else None
        if e is None:
            # zeta_M ** (z a) has order dividing rad(M); write it in zeta_N
            t = (z * a) % M
            g = gcd(t, M)
            order = M // g
            assert N % order == 0
            e = (t // g) * (N // order)
        if re:
            terms.append((re, e))
        if im:
            terms.append((im, e + N // 4))
    return K.from_exps(terms)


def delta(x, M, primes, K):
    """(P0, Q0): products of x_S over |S| even / odd."""
    P0, Q0 = K.one(), K.one()
    for r in range(len(primes) + 1):
        for S in itertools.combinations(primes, r):
            v = char_value(x, M, primes, S, K)
            if r % 2 == 0:
                P0 = K.mul(P0, v)
            else:
                Q0 = K.mul(Q0, v)
    return P0, Q0


# ---------------------------------------------------------------- 0. self-test
def selftest():
    assert cyclotomic(1) == [-1, 1]
    assert cyclotomic(4) == [1, 0, 1]
    assert cyclotomic(12) == [1, 0, -1, 0, 1]
    assert cyclotomic(9) == [1, 0, 0, 1, 0, 0, 1]
    K = Cyc(12)
    i = K.zeta(3)
    assert K.mul(i, i) == [-c for c in K.one()]
    assert K.norm(K.add(K.one(), i)) == 4  # N_{Q(zeta12)/Q}(1+i) = 2**2
    assert K.root_order(K.zeta(5)) == 12
    K3 = Cyc(3)
    assert K3.root_order([-c for c in K3.one()]) == 2
    assert K3.root_order(K3.zeta(1)) == 3
    # Delta(p**a X) = rho**(-1) for u = 3 (example in the note)
    x = {1: (3, 0)}
    P0, Q0 = delta(x, 9, [3], K3)
    d = K3.divide(P0, Q0)
    assert d is not None and K3.root_order(d) == 3
    assert gr_mul({1: (1, 0)}, {8: (1, 0)}, 9) == {0: (1, 0)}
    say("self-test passed")


# ---------------------------------------------------------------- A. identities
def part_A():
    rnd = random.Random(20261023)
    M = 15
    for trial in range(200):
        H = [{z: (rnd.choice((-1, 1)), 0) for z in range(M)} for _ in range(4)]
        def lin(s):
            return {z: (sum(s[j] * H[j][z][0] for j in range(4)), 0) for z in range(M)}
        # Klein characters (psi(a), psi(b)); h = H0 + a H1 + b H2 + ab H3
        c2 = lin((1, 1, 1, 1))
        d2 = [lin((1, -1, 1, -1)), lin((1, 1, -1, -1)), lin((1, -1, -1, 1))]
        for v in [c2] + d2:
            assert all(val[0] % 2 == 0 for val in v.values())
        bsum = {z: 0 for z in range(M)}
        for d in d2:
            for z in range(M):
                diff = d[z][0] - c2[z][0]          # = 2 (d - c) coefficientwise
                assert diff % 4 == 0               # d - c = -2 b with b integral
                bsum[z] += -diff // 4
        for z in range(M):
            assert bsum[z] == H[1][z][0] + H[2][z][0] + H[3][z][0]
            assert bsum[z] % 2 == 1               # b1+b2+b3 = J mod 2
        # cyclic 2-part: g = (H0-H2)/2 + i (H1-H3)/2 integral, and the identity
        # g - c = (i-1) b - (H2 + i H3), b = (H1+H3)/2
        for z in range(M):
            h0, h1, h2, h3 = (H[j][z][0] for j in range(4))
            assert (h0 - h2) % 2 == 0 and (h1 - h3) % 2 == 0 and (h1 + h3) % 2 == 0
            g = ((h0 - h2) // 2, (h1 - h3) // 2)
            c = (h0 + h1 + h2 + h3) // 2
            b = (h1 + h3) // 2
            rhs = (-b - h2, b - h3)              # (i-1) b - (h2 + i h3)
            assert (g[0] - c, g[1]) == rhs
            # h2 + i h3 = (1+i) + 2 U
            assert (h2 - 1) % 2 == 0 and (h3 - 1) % 2 == 0
    say("A: Klein and cyclic 2-part identities verified on 200 random sign arrays")


# ---------------------------------------------------------------- B. exhaustive
def vectors_sumsq(length, total, gaussian):
    """All vectors of the given length with sum of |x_z|**2 = total and entries
    in Z (real case) or Z[i] (Gaussian case)."""
    vals = {}
    for a in range(-3, 4):
        for b in (range(-3, 4) if gaussian else [0]):
            n = a * a + b * b
            if 0 < n <= total:
                vals.setdefault(n, []).append((a, b))
    out = []

    def rec(pos, left, cur):
        if pos == length:
            if left == 0:
                out.append(list(cur))
            return
        cur.append((0, 0))
        rec(pos + 1, left, cur)
        cur.pop()
        for n, lst in vals.items():
            if n <= left:
                for v in lst:
                    cur.append(v)
                    rec(pos + 1, left - n, cur)
                    cur.pop()
    rec(0, total, [])
    return out


def partial_vectors(length, total, gaussian):
    """Vectors of the given length with sum of |x_z|**2 <= total, grouped by norm."""
    by = {}
    for t in range(total + 1):
        for v in vectors_sumsq(length, t, gaussian):
            by.setdefault(t, []).append(v)
    return by


def solutions_norm9(gaussian):
    """All x in Z[C_9] (or Z[i][C_9]) with x x* = 9, by a split enumeration
    (5 + 4 coordinates); a float filter on the nine character values is used only
    to discard candidates, every survivor is then checked exactly."""
    import numpy as np
    M = 9
    w = np.exp(2j * np.pi * np.arange(M) / M)
    left = partial_vectors(5, 9, gaussian)
    right = partial_vectors(4, 9, gaussian)
    rarr = {}
    for t, lst in right.items():
        R = np.array([[complex(a, b) for a, b in v] for v in lst])
        # character values of the right part: positions 5..8
        E = np.array([[w[(k * (5 + j)) % M] for k in range(M)] for j in range(4)])
        rarr[t] = (lst, R @ E)
    EL = np.array([[w[(k * j) % M] for k in range(M)] for j in range(5)])
    sols = []
    ncand = 0
    for t, lst in left.items():
        if 9 - t not in rarr:
            continue
        rl, rv = rarr[9 - t]
        L = np.array([[complex(a, b) for a, b in v] for v in lst]) @ EL
        for i, v in enumerate(lst):
            tot = L[i][None, :] + rv
            ok = np.all(np.abs(np.abs(tot) ** 2 - 9.0) < 1e-6, axis=1)
            ncand += len(rl)
            for j in np.nonzero(ok)[0]:
                vec = v + rl[j]
                x = {z: vec[z] for z in range(M) if vec[z] != (0, 0)}
                if gr_mul(x, gr_star(x, M), M) == {0: (9, 0)}:
                    sols.append(x)
    return ncand, sols


def part_B():
    M = 9
    for gaussian in (False, True):
        K = Cyc(12) if gaussian else Cyc(3)
        ncand, sols = solutions_norm9(gaussian)
        orders = {}
        for x in sols:
            P0, Q0 = delta(x, M, [3], K)
            d = K.divide(P0, Q0)
            assert d is not None, "Delta not integral"
            o = K.root_order(d)
            assert o is not None and o in (1, 3, 9), ("bad order", o)
            # residue 1 above 3: (d - 1)/(1 - zeta_3) integral
            onem = K.sub(K.one(), K.zeta(K.N // 3))
            assert K.divide(K.sub(d, K.one()), onem) is not None
            orders[o] = orders.get(o, 0) + 1
        say("B: Z%s[C_9], x x* = 9: %d candidates, %d solutions, Delta orders %s"
            % ("[i]" if gaussian else "", ncand, len(sols), dict(sorted(orders.items()))))
        assert len(sols) > 0


# ---------------------------------------------------------------- C. necessity
D_KLEIN = [(0, 0, 0, 0), (0, 0, 0, 1), (0, 0, 0, 2), (0, 1, 0, 0), (0, 1, 0, 1),
           (0, 1, 1, 1), (0, 1, 1, 2), (0, 1, 2, 0), (0, 1, 2, 2), (1, 0, 0, 1),
           (1, 0, 1, 1), (1, 0, 2, 1), (1, 1, 0, 2), (1, 1, 1, 1), (1, 1, 2, 0)]
D_C4 = [(0, 0, 0, 0), (0, 0, 1, 0), (0, 0, 2, 0), (1, 0, 0, 1), (1, 0, 1, 2),
        (1, 0, 2, 0), (2, 0, 2, 0), (2, 0, 2, 1), (2, 0, 2, 2), (3, 0, 0, 1),
        (3, 0, 0, 2), (3, 0, 1, 0), (3, 0, 1, 1), (3, 0, 2, 0), (3, 0, 2, 2)]


def check_ds(D, mods):
    els = list(itertools.product(*[range(m) for m in mods]))
    cnt = {}
    for x in D:
        for y in D:
            if x != y:
                g = tuple((x[t] - y[t]) % mods[t] for t in range(len(mods)))
                cnt[g] = cnt.get(g, 0) + 1
    return len(D) == 15 and all(cnt.get(g, 0) == 6 for g in els if any(g))


def part_C():
    K = Cyc(3)
    assert check_ds(D_KLEIN, (2, 2, 3, 3)), "D_KLEIN is not a (36,15,6) difference set"
    assert check_ds(D_C4, (4, 1, 3, 3)), "D_C4 is not a (36,15,6) difference set"
    say("C: two (36,15,6) Hadamard difference sets verified (C2xC2xC3xC3 and C4xC3xC3)")
    # Klein: project C3xC3 -> C3 (first factor); the 3-part C3 has order 3 < 9 = u**2
    D = set(D_KLEIN)
    H = [{z: 0 for z in range(3)} for _ in range(4)]
    for a in range(2):
        for b in range(2):
            j = a + 2 * b
            for z in range(3):
                for w in range(3):
                    H[j][z] += -1 if (a, b, z, w) in D else 1
    def lin(s):
        return {z: (sum(s[j] * H[j][z] for j in range(4)) // 2, 0) for z in range(3)}
    c = lin((1, 1, 1, 1))
    ds = [lin((1, -1, 1, -1)), lin((1, 1, -1, -1)), lin((1, -1, -1, 1))]
    resid = []
    for x in [c] + ds:
        x = {k: v for k, v in x.items() if v != (0, 0)}
        assert gr_mul(x, gr_star(x, 3), 3) == {0: (9, 0)}
        P0, Q0 = delta(x, 3, [3], K)
        d = K.divide(P0, Q0)
        resid.append((d, K.root_order(d) if d is not None else None))
    say("C: Klein, projected to C2xC2xC3: Delta(c), Delta(d_i) =",
        [(r[0], r[1]) for r in resid])
    Rs = []
    for (d, o) in resid[1:]:
        Rs.append(K.divide(d, resid[0][0]) if d is not None and resid[0][0] is not None else None)
    say("C: R_i =", Rs, "(the argument needs all R_i = 1; it fails, as it must)")
    assert any(R != K.one() for R in Rs)


# ---------------------------------------------------------------- D. self-conjugacy
def part_D():
    K = Cyc(12)
    x = {0: (0, -1), 1: (1, 1), 2: (1, 1)}       # -i + (1+i) X + (1+i) X**2
    assert gr_mul(x, gr_star(x, 3), 3) == {0: (5, 0)}
    P0, Q0 = delta(x, 3, [3], K)
    d = K.divide(P0, Q0)
    say("D: x = -i + (1+i)X + (1+i)X**2 in Z[i][C_3], x x* = 5; x(1)/x(rho) integral?",
        d is not None)
    assert d is None
    # 5 is not self-conjugate modulo 12, but it is modulo 3
    assert not selfconj(5, 12) and selfconj(5, 3)


# ---------------------------------------------------------------- E. table
def selfconj(q, m):
    """True iff q**j = -1 (mod m) for some j (m >= 3), i.e. complex conjugation
    fixes the primes above q in Q(zeta_m); trivially true for m <= 2."""
    if m <= 2:
        return True
    t = q % m
    s = t
    for _ in range(m):
        if s == m - 1:
            return True
        s = s * t % m
        if s == t:
            break
    return False


def partitions(n, maxp=None):
    if maxp is None:
        maxp = n
    if n == 0:
        yield []
        return
    for k in range(min(n, maxp), 0, -1):
        for rest in partitions(n - k, k):
            yield [k] + rest


def part_E():
    total = cov_u = cov_c = lander_open = nolander = 0
    for u in range(3, 106, 2):
        f = factor(u)
        primes = sorted(f)
        choices = [list(partitions(2 * f[p])) for p in primes]
        for two in ("C4", "C2xC2"):
            for combo in itertools.product(*choices):
                total += 1
                cyc = [p for p, part in zip(primes, combo) if len(part) == 1]
                non = [p for p, part in zip(primes, combo) if len(part) > 1]
                if not cyc:
                    nolander += 1
                    continue
                # Corollary: some nonempty T within the cyclic primes such that every
                # other prime q | u is self-conjugate modulo rad(w_T) (4 rad(w_T))
                ok = False
                for k in range(1, len(cyc) + 1):
                    for T in itertools.combinations(cyc, k):
                        r = 1
                        for p in T:
                            r *= p
                        m = 4 * r if two == "C4" else r
                        if all(selfconj(q, m) for q in primes if q not in T):
                            ok = True
                if ok:
                    if non:
                        cov_c += 1
                    else:
                        cov_u += 1
                else:
                    lander_open += 1
    say("E: abelian groups of order 4u**2, u odd in [3,105]:", total)
    say("   no cyclic odd Sylow (Lander silent):", nolander)
    say("   covered, all odd Sylow cyclic (unconditional):", cov_u)
    say("   covered via self-conjugacy:", cov_c)
    say("   Lander hypothesis but not covered:", lander_open)
    assert (total, nolander, cov_u, cov_c, lander_open) == (424, 204, 104, 53, 63)
    return total



# ---------------------------------------------------------------- F. classical criteria
def lcm(a, b):
    return a // gcd(a, b) * b


def carmichael(fm):
    L = 1
    for p, e in fm.items():
        if p == 2:
            l = 1 if e == 1 else (2 if e == 2 else 2 ** (e - 2))
        else:
            l = (p - 1) * p ** (e - 1)
        L = lcm(L, l)
    return L


def mult_order(q, m, fm):
    if m == 1:
        return 1
    o = carmichael(fm)
    for p in factor(o):
        while o % p == 0 and pow(q, o // p, m) == 1:
            o //= p
    return o


def selfconj_fast(q, m, fm):
    if m <= 2:
        return True
    o = mult_order(q, m, fm)
    return o % 2 == 0 and pow(q, o // 2, m) == m - 1


def F_schmidt(m, n):
    """F(m,n) of Schmidt / Leung-Schmidt (field descent), Definition 2.6 of
    'The field descent method'."""
    fm = factor(m)
    D = sorted(factor(n))
    ps = sorted(fm)
    res = 1
    for p in ps:
        c = fm[p]
        for b in range(1, c + 1):
            ok = True
            for q in D:
                if q == p and (p, b) != (2, 1):
                    continue
                if b == c:
                    continue
                if q != p:
                    if m % 2 == 1 or q == 2:
                        mq = 1
                        for r in ps:
                            if r != q:
                                mq *= r
                    else:
                        mq = 4
                        for r in ps:
                            if r not in (2, q):
                                mq *= r
                    if pow(q, mult_order(q, mq, factor(mq)), p ** (b + 1)) != 1:
                        continue
                ok = False
                break
            if ok:
                break
        res *= p ** b
    return res


def phi(n):
    r = n
    for p in factor(n):
        r = r // p * (p - 1)
    return r


def part_F(UMAX=200000):
    """Which groups C2 x C_{2u**2} (u odd, not a prime power) survive two classical
    criteria: (i) the exponent bound under self-conjugacy (some p | u
    self-conjugate modulo the p-free part 2(u/p**a)**2 of the exponent excludes the
    group, since the cyclic Sylow p-subgroup has exponent p**(2a) > p**a);
    (ii) the field descent bound n <= |H| F**2/(4 phi(F)) for G = A x H, A cyclic,
    here A = C_{u**2}, H = C2 x C2, F = gcd(u**2, F(2u**2, u**2)).
    For comparison, the circulant bound u phi(u) <= F(u**2, u) is also evaluated."""
    for m in range(3, 300):
        for q in (3, 5, 7, 11, 13, 17, 19, 23):
            if gcd(q, m) == 1:
                assert selfconj(q, m) == selfconj_fast(q, m, factor(m))
    klein, circ = [], []
    for u in range(3, UMAX + 1, 2):
        f = factor(u)
        if len(f) == 1:
            continue
        excl_k = excl_c = False
        for p in f:
            for two, store in ((1, "k"), (2, "c")):
                fm = {2: two}
                for q in f:
                    if q != p:
                        fm[q] = 2 * f[q]
                m = 1
                for a_, e in fm.items():
                    m *= a_ ** e
                if selfconj_fast(p, m, fm):
                    if store == "k":
                        excl_k = True
                    else:
                        excl_c = True
        if not excl_c and u * phi(u) <= F_schmidt(u * u, u):
            circ.append(u)
        if not excl_k:
            Fv = gcd(u * u, F_schmidt(2 * u * u, u * u))
            if u * u * phi(Fv) <= Fv * Fv:
                klein.append(u)
    say("F: u odd, not a prime power, u <= %d, not excluded by the exponent bound"
        " and the field descent bound:" % UMAX)
    say("   circulant C_{4u**2}:   ", circ)
    say("   Klein C2 x C_{2u**2}:  ", klein)
    assert circ and circ[0] == 11715  # agrees with Leung-Schmidt, v < 548964900
    assert klein == circ == [11715, 16401, 82005]
    return klein


if __name__ == "__main__":
    selftest()
    part_A()
    part_B()
    part_C()
    part_D()
    part_E()
    part_F()
    say("all checks passed in %.1f s" % (time.time() - T0))