"""dualrs.py -- the chain law's depth side through a class the chain's
step keeps: does a hidden window always hold a k with q g_k = 0?

QUESTION. dualn.py's setting: f monic of degree n, f(0) = +/-q, O =
Z[x]/(f), G = O/I, g_k = beta^k mod I, D = {0, ..., q - 1}, T_W the
sums over k in W of D g_k. The class RS: a finite abelian G with a
sequence g_0, g_1, ... generating it and satisfying
sum_m a_m g_(k+m) = 0 for every k >= 0 (a_m f's coefficients) and
g_k = 0 for k >= some t. O/I read at depth t is in RS, and so is
G/<g_k> with the images. GC: in RS, if T_W is uniform on G != 0, then
q g_k = 0 for some k in W.

HAND (frozen before the runs).
  H1 If q g_k = 0, D g_k is q/r copies of <g_k>, r its order, so T_W is
     uniform on G iff T_(W - k) is on G/<g_k>, and a chain step in the
     quotient is one in G; GC therefore gives the depth side by
     induction on |G| + |W|. Uniformity is monotone in W, so GC iff
     W* = {k : q g_k != 0} never hides G != 0.
  R1 q = l^e. U = O/beta^t O is an l-group, so U = O0/beta^t O0 with
     O0 = Z_l[x]/(f0), f0 the factor of f over Z_l with f0 = x^r mod l
     (Hensel; beta is a unit on the other factor). So a G in RS(f)
     carries a sequence obeying f0's recurrence, f0 = x^r + a_(r-1)
     x^(r-1) + ... + a_0, l | a_i, v_l(a_0) = e, read mod exp(G); and
     every such f0 mod exp(G) is the factor, mod exp(G), of an integer
     f with f(0) = +/-q, f = f0 (x + c) with c a unit and the constant
     term set to +/-q. One sweep over f0 of degree r covers every
     integer base whose local factor has degree r, and only those. In
     an l-group each f0-sequence reaches 0 (g_(k+r) lies in
     l <g_k, ..., g_(k+r-1)>), and G has rank <= r.
  R2 A cyclic G never hides: a faithful chi has q chi(g_k) != 0 at
     every k in W*. Skipped, with every G of exponent <= q (W* empty).
  R3 g_k -> v^k g_k, v a unit of Z_l, keeps every order, so W* and
     hiding, and maps f0 to a_i v^(r-i). At l = 2, r = 3 cubing is onto
     the units: a_0 = q. At l = 3, r = 3: a_0 = 9u, u in {1, 4, 7}.
  R4 chi in B_k iff chi in (q g_k)^perp - g_k^perp; W* hides iff the
     B_k, q g_k != 0, cover every chi != 0. One bitmask per element.
  R7 On (Z/l^m)^d, GL_d(Z/l^m) is transitive on the elements of each
     order and commutes with the recurrence, so g_0 runs over 0 and
     (l^j, 0, ..., 0) only.

PREDICTIONS (frozen). Every initial r-tuple, run until it reaches 0.
  K  integer cubics f = x^3 + a2 x^2 + a1 x + c0, c0 = +/-q, a1, a2 in
     -2..2, G every product of cyclic l-groups of order <= 16 (q = 2,
     4) or 27 (q = 3), and q = 4 at order 32 with a1, a2 even, every
     triple of G^3 run to 64 terms, kept iff it reaches 0: 0 hiding.
  A  q = 4, r = 3, G of order 64, rank <= 3, exp >= 8, every a_1, a_2
     even mod exp(G), a_0 = 4;
  B  q = 8, r = 3, G of order 64, exp >= 16, a_1, a_2 even, a_0 = 8;
  C  q = 9, r = 3, G of order 81, exp >= 27, a_1, a_2 in 3Z mod exp,
     a_0 = 9u, u in {1, 4, 7};
  D  q = 4, r = 4, G of order <= 32, exp >= 8, a_1..a_3 even mod exp,
     a_0 = 4u, u odd;
  E  the groups every invariant factor of which exceeds q, past order
     64: q = 4 at (16, 8), (16, 16), (8, 8, 8), (32, 32); q = 8 at
     (16, 16); q = 9 at (27, 27), a_0 = 9.
  Predicted 0 hiding everywhere. KILL: any hiding sequence, printed
  with G, f0 and its r-tuple; GC then dies.
  L2 (frozen with T2 below, before its code) every noncyclic rank-2
     l-group of order <= 256 at q = 2, 4, 8 and <= 243 at q = 3, 9,
     f0 = x^2 + a_1 x + q u (u over the units mod squares, a_1 in lZ
     mod exp G), and (added after the first run, whose q = 9 arm held
     no group with both factors above q) (32, 32) at q = 4, 8 and
     (27, 27) at q = 9; every (g_0, g_1) run to 0: at most 2 distinct lines
     <q g_k> with ord q g_k = l. KILL: a sequence with 3. Its own
     print: the sequences at 2, which must be > 0 for the bound to be
     tested.
  S2 r = 2 with no nilpotence: q = 2, 4 at order 16, 32, 64, q = 3 at
     27, 81, every group, every a_0 = q u and every a_1 mod exp G, W*
     cut to [0, 8) and [0, 12): 0 hiding (T2 reads any window).
CONTROLS (read first).
  C1 "reaches 0" dropped, W* cut to [0, 8), q = 2, order 16, integer
     cubics: hiding sequences exist.
  C2 O/I at every ideal read at some depth, cubic bases of dualn.py at
     q^J <= 128, W* over [0, t): no hiding.
  C3 R7's cut at (4, 4), C1's only hiding group: one g_0 per order,
     and each hides iff its order class does.
  CB the nil arm, which no hiding reaches: the engine's witness count
     (characters no B_k covers) against a brute-force count over every
     character and every k in W*, at sampled eventually-zero sequences
     of (8, 8), (16, 8) at q = 4 and (27, 9) at q = 9, and at r = 2 of
     (8, 8) at q = 4 and (9, 9) at q = 3; equal at every one.
  C4 L2's code at r = 3: (8, 8), q = 4, x^3 + 4, every triple: some
     sequence puts q g_k on all 3 lines (the example closing DERIVED).
     Below 3, L2 is void.
  C5 the line counter itself: q = 1, f = x^7 on (Z/2)^3, seven free
     terms: some sequence hits all 7 lines (a counter reading only the
     first l + 1 line bits stops at 3).

FINDINGS (entered after the run, from its printed output).
  C1 6144 hiding sequences without nilpotence. C2 0 of 1948 ideals
     read at some depth hidden by their W*. C3 12 hiding (pair, order)
     cells, 0 off, every order held. CB 1080 sequences, 0 off, 31
     distinct witness counts, the least 12. C4 24576 triples of (8, 8)
     hit all 3 lines. C5 most 7 lines.
  L2 most 2 lines, 3095568 sequences at 2 under R7's cut: q 2, 861
     (G, f0) pairs, 2047632 at 2; q 4, 849, 640640; q 8, 689, 13312;
     q 3, 103, 385236; q 9, 115, 8748. S2 26392 runs, 0 hiding.
  K  0 hiding: kept 419112 (q 2, order <= 16), 144012 (q 3, <= 27),
     419112 (q 4, <= 16), 4128768 (q 4, order 32, a1, a2 even); each
     count holds every initial triple reaching 0, the zero sequence,
     the non-generating ones (which cannot hide) and automorphic
     copies among them.
  A  0 hiding: 416 (G, f0) pairs, 105119744 sequences under R7's cut.
  B  0 hiding: 384 pairs, 100663296 sequences.
  C  0 hiding: 81 pairs, 43046721 sequences (at exp 27, a_0 = 9u reads
     u mod 3 only, so the three u are one).
  D  (--full) 0 hiding: 1216 pairs, 1212153856 sequences.
  E  (--full) 0 hiding: q = 4, 400 pairs, 1782579200 sequences under
     R7's cut; q = 8, 64 pairs, 20971520; q = 9, 81 pairs, 172186884.
  The kill misses at every family: GC holds at this scale, a pattern
  at r >= 3; at r <= 2 it is T2 below.

DERIVED (by hand: GC at r <= 2, T2; toward r >= 3).
  D1 RS is the quotients of U = O/beta^t O by its subgroups H, g_k the
     class of beta^k: U is generated by 1, ..., beta^(t-1) under f's
     relations and beta^t = 0. T_t is a complete residue system of U
     (O/beta O = Z/q), so with Z = [0, t) - W*, T_(W*) + T_Z = U once
     each, and W* hides U/H iff T_(W*) + H covers U m = |H|/q^|Z| times.
  D2 V_k = <g_k, g_(k+1), ...> has V_k/V_(k+1) cyclic on g_k of order
     n_k | q, as in the ring. A character built down from V_(t-1) meets
     a constraint only at k in W* with n_k < q: there q chi(g_k) =
     (q/n_k) chi(n_k g_k) is fixed before chi(g_k) is chosen.
  D3 The filtration alone does not carry GC: a no-depth miss's g_0,
     ..., g_5 with elements z, q z = 0, appended to span the q g_k
     satisfies q y_i in <y_(i+1), ...> for every i, and its window
     still hides. A proof must use f's recurrence past u = max W*.
  D4 In a least counterexample every nonzero subgroup K has a chain
     over W* started at K that reaches G (the quotient G/K is smaller),
     so every subgroup of order l is <q g_k> for some k in W*.
  R5 S <= G is CLOSED when k in W*, q g_k in S imply g_k in S. A
     nonzero proper closed S reduces GC at G to GC at G/S (in RS): a
     witness chi of G/S for W*(G/S) lifts, since at k in W* with
     q g_k in S, g_k is in S and chi(g_k) = 0. So a least
     counterexample has none (D4 is |S| = l).
  R6 (case A) If an invariant factor of G has order <= q, (qG)[l] is
     not G[l], and a line of G[l] outside qG is closed (no q g_k lies
     in it). So a least counterexample is CASE B: every invariant
     factor exceeds q, every line of G[l] is <q g_k> with o_k = e + 1
     (ord g_k = l^(o_k)).
  R8 The closure of qG (the least closed subgroup over it) is qG +
     <g_k : k in W*>: every k in W* enters at once. So in a least
     counterexample the g_k, k in W*, generate G.
  R9 T_t = T_(W*) + T_Z is a complete residue system of U (D1), so in
     G, mu_(W*) * mu_Z = |H| u_G (mu_S the image of T_S counted with
     multiplicity, u_G uniform). For k in Z, D g_k covers <g_k>
     q/ord(g_k) times, and a sum of uniform measures on subgroups is
     uniform on their sum: mu_Z is a multiple of u_(G_Z), G_Z =
     <g_k : q g_k = 0>. So W* hides G iff mu_(W*) is invariant under
     translation by G_Z, and a witness is nonzero on G_Z.
  R10 (the dual) U^ is the f0-sequences c in Q_l/Z_l with c_k = 0 for
     k >= t, G^ a subgroup of it, c_k = chi(g_k). q c_k is fixed by
     c_(k+1), ..., c_(k+r) (a_0 = q times a unit), so the last nonzero
     term of every c is killed by q, and a witness's last nonzero term
     sits off W*.
  R11 (r = 2) For a lattice M in K = O0 (x) Q_l, the beta^k (k >= 0) in
     M - lM put their images on at most two lines of M/lM. Write
     beta^m = A_m + B_m beta: A_1 = 0 and A_(m+1) = -a_0 B_m, so l | A_m
     for m >= 1. Let i be the least such k. If some such j has its
     image off i's line, M = Z_l beta^i + Z_l beta^j (Nakayama), a
     Q_l-basis of K (M is a full lattice), so B_d != 0, d = j - i; any
     other such k, m = k - i, has beta^k = x beta^i + y beta^j with
     y = B_m/B_d and x = A_m - y A_d in Z_l and l | x, so it lies on j's
     line. L2 prints the bound; at r = 3 beta^m mod l ranges over
     <beta, beta^2>, and G = Z/8 x Z/8, f0 = x^3 + 4, g = (1,0), (0,1),
     (1,1) puts q g_0, q g_1, q g_2 on all three lines (C4).
  T2 (GC at r <= 2) q = l^e, f0 monic of degree <= 2 over Z_l,
     l | f0(0), G a finite l-group, g_k obeying f0's recurrence, W
     finite: W* never hides G != 0; no nilpotence is used (S2 sweeps
     it). If the g_k span a proper subgroup, a chi != 0 vanishing on
     it is a witness; else G = O0/L, L the kernel of x^k -> g_k (a
     Z_l-lattice, not always an ideal), rank <= 2. R5 and R6 read
     only G, the g_k and W*, and a quotient keeps f0 and W, so a
     least counterexample is case B of rank 2 (R2), each of the l + 1
     >= 3 lines of G[l] some <q g_k>. Both invariant factors above q
     put L in l^(e+1) O0; with M = l^-(e+1) L, x -> l^e x maps M/lM
     onto G[l] = l^e M / l^(e+1) M, and <q g_k> is a line iff beta^k
     lies in M - lM, the line of beta^k mod lM. R11 allows two. No
     invertibility of lattices enters, so a double root of f0 costs
     nothing. With H1, the depth side holds wherever f's local factor
     has degree <= 2.
  Case B stays open at r >= 3. The aim, seen and not printed here: at
  r = 3, over the cubic local factors (a_0 = q, a_1, a_2 even), every
  generating eventually-zero sequence of (8, 8), (16, 8), (16, 16),
  (8, 8, 8) at q = 4 and (16, 16) at q = 8, under R7's cut, had a
  witness chi with chi(g_s) != 0, s the sequence's last k with
  g_k != 0, and one of order exp G, each found apart. The sequences
  with every line of G[l] some <q g_k> have q^|W*| = |G| at (8, 8),
  where Hajos's theorem kills a hiding, and q^|W*| > |G| at (16, 8) at
  q = 4 and (16, 16) at q = 4, 8, where it does not reach.

RUN RECORD. Default 12/12, 86.4 s, 68.3 MB; D 288.3 s, 67.6 MB, under a
memory guard. E last ran before the sweep was split into block(), the
split leaving A, B, C and K's counts unchanged: 1797.8 s, peak 218.3
MB, nearly all of it (Z/32)^2.
"""

import os
import random
import sys
import time
from itertools import product

os.environ.setdefault("OPENBLAS_NUM_THREADS", "1")
import numpy as np  # noqa: E402

CHECKS = []
CHUNK = 1 << 18


def check(name, ok, detail=""):
    CHECKS.append(bool(ok))
    print(f"  {'PASS' if ok else 'FAIL'}  {name}" +
          (f"  ({detail})" if detail else ""), flush=True)


def groups(l, top, exact=False):
    """Every product of cyclic l-groups of order <= top (== top if
    exact), as moduli lists, largest first."""
    out = []

    def rec(rest, cap, cur):
        if cur and (not exact or rest == 0):
            out.append(tuple(l ** x for x in cur))
        for e in range(min(cap, rest), 0, -1):
            rec(rest - e, e, cur + [e])
    e = 0
    while l ** (e + 1) <= top:
        e += 1
    rec(e, e, [])
    return out


class Group:
    """G = prod Z/m_i, elements as indices 0..N-1 (0 the identity)."""

    def __init__(self, ms):
        self.ms = ms
        self.E = max(ms)
        self.els = np.array(list(product(*[range(m) for m in ms])),
                            dtype=np.int64)
        self.N = N = len(self.els)
        st = np.cumprod([1] + [m for m in ms[::-1]])[:-1][::-1]
        self.st = np.array(st, dtype=np.int64)
        self.add = self.index(self.els[:, None, :] +
                              self.els[None, :, :]).ravel()
        R = self.els * np.array([self.E // m for m in ms], dtype=np.int64)
        zero = (self.els @ R.T) % self.E == 0          # [x, chi]
        self.W = (N + 63) // 64
        bits = np.zeros((N, self.W * 64), dtype=bool)
        bits[:, :N] = zero
        self.perp = np.packbits(bits.reshape(N, self.W, 64)[:, :, ::-1],
                                axis=2, bitorder="big") \
            .view(">u8").astype(np.uint64).reshape(N, self.W)
        full = np.ones(self.W * 64, dtype=bool)
        full[N:] = False
        full[0] = False                                 # chi = 0
        self.full = np.packbits(full.reshape(self.W, 64)[:, ::-1],
                                axis=1, bitorder="big") \
            .view(">u8").astype(np.uint64).reshape(self.W)

    def index(self, v):
        return ((v % np.array(self.ms, dtype=np.int64)) @ self.st)

    def mul(self, c):
        return self.index(c * self.els)


def homocyclic_firsts(G):
    """R7: on (Z/l^m)^d, g_0 in {0, (l^j, 0, ..., 0)} up to Aut."""
    m, d = G.ms[0], len(G.ms)
    if any(x != m for x in G.ms):
        return None
    out, c = [0], 1
    while c < m:
        out.append(int(G.index(np.array([c] + [0] * (d - 1)))))
        c *= [p for p in range(2, m + 1) if m % p == 0][0]
    return out


def block(G, q, a, ids, nil=True, window=None, terms=64):
    """Run the sequences whose initial r-tuples are the base-N digits
    of ids, g_(k+r) = -(a_0 g_k + ... + a_(r-1) g_(k+r-1)). Returns
    (cover, keep): cover the B_k's union as bitmasks over chi, keep
    whether the sequence reaches 0 (the last r of `terms` zero; run
    until all reach 0 when nil); window: cover only k < window."""
    r, N = len(a), G.N
    mula = [G.mul(-c) for c in a]
    qm = G.mul(q)
    seq = [ids // N ** (r - 1 - j) % N for j in range(r)]
    cover = np.zeros((len(ids), G.W), dtype=np.uint64)
    k = 0
    while True:
        g = seq[k]
        if window is None or k < window:
            qg = qm[g]
            m = G.perp[qg] & ~G.perp[g]
            m[qg == 0] = 0
            cover |= m
        nxt = mula[0][seq[k]]
        for j in range(1, r):
            nxt = G.add[nxt * N + mula[j][seq[k + j]]]
        seq.append(nxt)
        k += 1
        if nil and window is None and k >= r and \
                not any(np.any(s) for s in seq[k:k + r]):
            break
        if k + r >= terms:
            break
    if nil:
        keep = np.all(np.stack([seq[j] == 0 for j in
                                range(len(seq) - r, len(seq))]), axis=0)
    else:
        keep = np.ones(len(ids), dtype=bool)
    return cover, keep


def sweep(G, q, a, nil=True, window=None, terms=64, firsts=None):
    """Every r-tuple of G (g_0 in firsts if given). Returns (kept,
    hiding, first hiding tuple)."""
    r, N = len(a), G.N
    rest = N ** (r - 1)
    blocks = [(f * rest, (f + 1) * rest) for f in firsts] if firsts \
        else [(0, N ** r)]
    kept = hid = 0
    first = None
    for b0, hi in blocks:
        for lo in range(b0, hi, CHUNK):
            ids = np.arange(lo, min(hi, lo + CHUNK), dtype=np.int64)
            cover, keep = block(G, q, a, ids, nil, window, terms)
            hides = np.all(cover == G.full, axis=1) & keep
            kept += int(keep.sum())
            hid += int(hides.sum())
            if first is None and hides.any():
                i = int(ids[np.argmax(hides)])
                first = tuple(tuple(G.els[i // N ** (r - 1 - j) % N])
                              for j in range(r))
    return kept, hid, first


def local_coeffs(l, E, n):
    """Every (a_1..a_n) with l | a_i, read mod E."""
    return list(product(range(0, E, l), repeat=n))


def integer_cubics(q, even=False):
    rng = (-2, 0, 2) if even else range(-2, 3)
    for c0 in (q, -q):
        for a1 in rng:
            for a2 in rng:
                yield (c0, a1, a2)


def popcount(x):
    """Bits set per row of a uint64 array."""
    v = x.view(np.uint8)
    return np.unpackbits(v, axis=-1).reshape(x.shape[0], -1).sum(axis=1)


def brute_witnesses(ms, q, a, tup):
    """Characters chi != 0 with chi(g_k) = 0 or q chi(g_k) != 0 at every
    k with q g_k != 0, by direct arithmetic; the sequence run to 0."""
    E, r = max(ms), len(a)
    s = [list(x) for x in tup]
    while any(any(x) for x in s[-r:]):
        s.append([(-sum(a[j] * s[-r + j][i] for j in range(r))) % m
                  for i, m in enumerate(ms)])
    W = [x for x in s if any((q * xi) % m for xi, m in zip(x, ms))]
    n = 0
    for c in product(*[range(m) for m in ms]):
        if not any(c):
            continue
        vals = [sum(ci * xi * (E // m) for ci, xi, m in zip(c, x, ms)) % E
                for x in W]
        n += all(v == 0 or (q * v) % E for v in vals)
    return n


def control():
    print("C  controls", flush=True)
    t0 = time.time()
    hid = 0
    for ms in groups(2, 16, exact=True):
        G = Group(ms)
        for a in integer_cubics(2):
            hid += sweep(G, 2, a, nil=False, window=8, terms=11)[1]
    check("C1 no nilpotence, W* cut to [0, 8), q 2 order 16: hiding "
          "exists", hid > 0, f"{hid}, {time.time() - t0:.1f} s")
    G = Group((4, 4))
    order = [int(np.lcm.reduce([m // np.gcd(int(x), m) for x, m in
                                zip(e, G.ms)])) for e in G.els]
    cells = off = 0
    for a in integer_cubics(2):
        for f in homocyclic_firsts(G):
            orbit = [x for x in range(G.N) if order[x] == order[f]]
            whole = sweep(G, 2, a, nil=False, window=8, terms=11,
                          firsts=orbit)[1]
            rep = sweep(G, 2, a, nil=False, window=8, terms=11,
                        firsts=[f])[1]
            cells += whole > 0
            off += (whole > 0) != (rep > 0)
    span = {order[f] for f in homocyclic_firsts(G)} == set(order)
    check("C3 R7: one representative per order of g_0, and it hides iff "
          "its orbit does, (4, 4)", span and cells > 0 and off == 0,
          f"{cells} hiding (pair, order) cells, {off} off, every order "
          f"{span}")
    rnd = random.Random(7)
    pairs = [((8, 8), 4, (4, 0, 0)), ((8, 8), 4, (4, 2, 6)),
             ((16, 8), 4, (4, 6, 8)), ((16, 8), 4, (12, 2, 4)),
             ((27, 9), 9, (9, 3, 6)), ((8, 8), 4, (4, 2)),
             ((8, 8), 4, (12, 6)), ((9, 9), 3, (3, 3)), ((9, 9), 3, (6, 0))]
    seen = bad = 0
    counts = set()
    for ms, q, a in pairs:
        G = Group(ms)
        r = len(a)
        ids = np.array(rnd.sample(range(G.N ** r), 3000), dtype=np.int64)
        cover, keep = block(G, q, a, ids, terms=10 ** 6)
        eng = popcount(G.full & ~cover)
        for i in np.where(keep)[0][:120]:
            tup = [tuple(int(x) for x in G.els[int(ids[i]) //
                                              G.N ** (r - 1 - j) % G.N])
                   for j in range(r)]
            b = brute_witnesses(ms, q, a, tup)
            seen += 1
            bad += b != int(eng[i])
            counts.add(b)
    check("CB the nil arm's witness count matches brute force",
          seen > 0 and bad == 0 and len(counts) > 3,
          f"{seen} sequences, {bad} off, {len(counts)} distinct counts, "
          f"least {min(counts)}, {time.time() - t0:.1f} s")
    return time.time() - t0


def section_c2():
    from dualdeg import powers
    from dualn import Base, cubic_bases
    from dualmult import depth
    t0 = time.time()
    n_id = hid = 0
    for a in cubic_bases():
        B = Base(a, 128)
        q = B.q
        for L in B.Is:
            t = depth(L, a, B.pw[1])
            if t is None or L.N == 1:
                continue
            n_id += 1
            pw = powers(a, t + 1)
            gk = [pw[k] for k in range(t)]
            Ws = [k for k in range(t)
                  if any(L.key(tuple(q * x for x in gk[k])))]
            chars = L.characters()[1:]
            N = L.N
            ok = True
            for c in chars:
                vals = [sum(ci * xi for ci, xi in zip(c, gk[k])) % N
                        for k in Ws]
                if not any(v and (q * v) % N == 0 for v in vals):
                    ok = False
                    break
            hid += ok
    check("C2 no ring ideal read at some depth is hidden by its W*",
          n_id > 0 and hid == 0,
          f"{hid} of {n_id}, {time.time() - t0:.1f} s")


def section_k():
    print("K  integer cubics, every triple to 64 terms", flush=True)
    t0 = time.time()
    total = 0
    for q, l, top, even, exact in [(2, 2, 16, False, False),
                                   (3, 3, 27, False, False),
                                   (4, 2, 16, False, False),
                                   (4, 2, 32, True, True)]:
        kept = hid = 0
        for ms in groups(l, top, exact=exact):
            G = Group(ms)
            for a in integer_cubics(q, even):
                k_, h_, _ = sweep(G, q, a)
                kept += k_
                hid += h_
        total += hid
        print(f"    q {q}, order {'' if exact else '<= '}{top}"
              f"{', a1, a2 even' if even else ''}: kept {kept}, hiding "
              f"{hid}, {time.time() - t0:.1f} s", flush=True)
    check("K  no eventually-zero sequence has W* hiding G", total == 0,
          f"{total}")


def line_bits(G, l):
    """Element index -> the bit of its line of G[l] (0 off G[l] - 0),
    and the number of lines."""
    lm, muls = G.mul(l), [G.mul(c) for c in range(1, l)]
    bit = np.zeros(G.N, dtype=np.int64)
    n = 0
    for x in range(1, G.N):
        if lm[x] == 0 and bit[x] == 0:
            for m in muls:
                bit[m[x]] = 1 << n
            n += 1
    return bit, n


def lines_hit(G, q, a, l, firsts=None):
    """Every r-tuple (g_0 in firsts if given), run to 0: per sequence
    the number of distinct lines <q g_k> with ord q g_k = l, as a
    histogram {count: sequences}."""
    r, N = len(a), G.N
    bit, n = line_bits(G, l)
    bit = bit[G.mul(q)]
    mula = [G.mul(-c) for c in a]
    rest = N ** (r - 1)
    blocks = [(f * rest, (f + 1) * rest) for f in firsts] if firsts \
        else [(0, N ** r)]
    hist = {}
    for b0, hi in blocks:
        for lo in range(b0, hi, CHUNK):
            ids = np.arange(lo, min(hi, lo + CHUNK), dtype=np.int64)
            seq = [ids // N ** (r - 1 - j) % N for j in range(r)]
            acc = np.zeros(len(ids), dtype=np.int64)
            k = 0
            while any(np.any(x) for x in seq[k:k + r]):
                acc |= bit[seq[k]]
                nxt = mula[0][seq[k]]
                for j in range(1, r):
                    nxt = G.add[nxt * N + mula[j][seq[k + j]]]
                seq.append(nxt)
                k += 1
            cnt = sum((acc >> b) & 1 for b in range(n))
            for c, m in zip(*np.unique(cnt, return_counts=True)):
                hist[int(c)] = hist.get(int(c), 0) + int(m)
    return hist


def section_l2():
    print("L2 lines <q g_k> of order l, r = 2", flush=True)
    t0 = time.time()
    h7 = lines_hit(Group((2, 2, 2)), 1, (0,) * 7, 2, firsts=[1])
    check("C5 the counter reads all 7 lines of (2, 2, 2): q = 1, f = x^7, "
          "seven free terms", max(h7) == 7, f"most {max(h7)}")
    h3 = lines_hit(Group((8, 8)), 4, (4, 0, 0), 2)
    check("C4 r = 3, (8, 8), q 4, x^3 + 4: some sequence hits all 3 "
          "lines", max(h3) == 3, f"histogram {h3}")
    top = two = 0
    for l, q, order, us, extra in [(2, 2, 256, (1, 3, 5, 7), []),
                                   (2, 4, 256, (1, 3, 5, 7), [(32, 32)]),
                                   (2, 8, 256, (1, 3, 5, 7), [(32, 32)]),
                                   (3, 3, 243, (1, 2), []),
                                   (3, 9, 243, (1, 2), [(27, 27)])]:
        pairs, hist = 0, {}
        for ms in [x for x in groups(l, order) if len(x) == 2] + extra:
            G = Group(ms)
            firsts = homocyclic_firsts(G)
            for a0 in sorted({q * u % G.E for u in us}):
                for (a1,) in local_coeffs(l, G.E, 1):
                    for c, n in lines_hit(G, q, (a0, a1), l,
                                          firsts).items():
                        hist[c] = hist.get(c, 0) + n
                    pairs += 1
        top = max(top, max(hist))
        two += hist.get(2, 0)
        print(f"    q {q}, rank 2, order <= {order}"
              f"{''.join(f', {x}' for x in extra)}: {pairs} (G, f0) pairs, "
              f"sequences by lines hit {dict(sorted(hist.items()))}, "
              f"{time.time() - t0:.1f} s", flush=True)
    check("L2 r = 2: no sequence puts q g_k on 3 lines", top <= 2 and
          two > 0, f"most {top}, {two} sequences at 2 under R7's cut")


def section_s2():
    t0 = time.time()
    scan = runs = 0
    for q, l, tops in [(2, 2, (16, 32, 64)), (4, 2, (16, 32, 64)),
                       (3, 3, (27, 81))]:
        for top in tops:
            for ms in groups(l, top, exact=True):
                G = Group(ms)
                for u in range(1, G.E):
                    if u % l == 0:
                        continue
                    for a1 in range(G.E):
                        for win in (8, 12):
                            scan += sweep(G, q, (q * u % G.E, a1), nil=False,
                                          window=win, terms=win + 3)[1]
                            runs += 1
    check("S2 r = 2 without nilpotence, every a_1, W* cut to [0, 8) and "
          "[0, 12): no hiding", runs > 0 and scan == 0,
          f"{runs} runs, {scan} hiding, {time.time() - t0:.1f} s")


def family(name, l, q, r, orders, a0s, only=None):
    t0 = time.time()
    seqs = hid = pairs = 0
    for top in orders:
        for ms in only or groups(l, top, exact=True):
            if len(ms) == 1 or len(ms) > r or max(ms) <= q:
                continue
            G = Group(ms)
            firsts = homocyclic_firsts(G)
            for a0 in sorted({a % G.E for a in a0s}):
                for rest in local_coeffs(l, G.E, r - 1):
                    kept, h, first = sweep(G, q, (a0,) + rest,
                                           terms=10 ** 6, firsts=firsts)
                    seqs += kept
                    pairs += 1
                    if h:
                        hid += h
                        print(f"    HIDES  G {ms}, f0 coeffs "
                              f"{(a0,) + rest}, {h} tuples, first {first}",
                              flush=True)
    print(f"  {name}: {pairs} (G, f0) pairs, {seqs} sequences, hiding "
          f"{hid}, {time.time() - t0:.1f} s", flush=True)
    check(f"{name.split()[0]}  no hiding", hid == 0, f"{hid}")


def main():
    t0 = time.time()
    full = "--full" in sys.argv
    control()
    section_c2()
    section_l2()
    section_s2()
    section_k()
    family("A q 4 cubic, order 64", 2, 4, 3, [64], [4])
    family("B q 8 cubic, order 64", 2, 8, 3, [64], [8])
    family("C q 9 cubic, order 81", 3, 9, 3, [81], [9, 36, 63])
    if full:
        family("D q 4 quartic, order <= 32", 2, 4, 4, [16, 32],
               [4 * u for u in range(1, 16, 2)])
        family("E q 4, every invariant factor above q", 2, 4, 3, [0], [4],
               only=[(16, 8), (16, 16), (8, 8, 8), (32, 32)])
        family("E q 8, (16, 16)", 2, 8, 3, [0], [8], only=[(16, 16)])
        family("E q 9, (27, 27)", 3, 9, 3, [0], [9], only=[(27, 27)])
    print(f"{sum(CHECKS)}/{len(CHECKS)} checks pass, "
          f"{time.time() - t0:.1f} s")
    return 0 if all(CHECKS) else 1


if __name__ == "__main__":
    sys.exit(main())