"""dualpw.py -- the primitive witness at q = 2^e: at every cubic local
factor, does every rank-2 2-group of type (2^N, 2^M), M > e, hold a
closed <h> at every generating sequence?

QUESTION. dualrs.py's setting at l = 2, q = 2^e, r = 3: f0 = x^3 + a2
x^2 + a1 x + a0 over Z_2 with a1, a2 even and v(a0) = e (R1), scaled to
a0 = 2^e (R3: cubing is onto the units of Z_2, and beta -> w beta keeps
a1, a2 even). G = Z/2^N x Z/2^M, e + 1 <= M <= N, g_k obeying f0's
recurrence mod 2^N and generating G. PW (dualrs.py R15): some <h>, h
not in 2G, is closed: no k with q g_k != 0 has q g_k in <h> and g_k
outside it. By R2 and T3 a least counterexample to GC at a cubic factor
is rank 2, by R14 both invariant factors above q and closure-free; a
closed <h> refutes that, and a quotient's exponent is at most G's, so
PW at every such G of exponent at most 2^N gives GC at q = 2^e, at
every cubic factor, on every G of exponent at most 2^N.

HAND (properties). O0 = Z_2[beta], beta^3 = -(a2 beta^2 + a1 beta +
2^e), O0/2^N O0 coordinates in 1, beta, beta^2, c_2 the beta^2
coordinate. At e = 1 f0 is Eisenstein and beta a uniformizer; at e = 2
it is not (v(a1) = 1: a linear times an Eisenstein quadratic factor;
v(a1) >= 2: irreducible, Z_2[beta] not maximal), at e = 3 it can split
(x^3 + 8 = (x + 2)(x^2 - 2x + 4)), and the engine reads
no ring property of O0: the pairing c_2(y x) has Gram matrix [[0, 0,
1], [0, 1, -a2], [1, -a2, a2^2 - a1]] on 1, beta, beta^2, det -1, so it
is perfect at every monic cubic, and x^k -> g_k is Z_2-linear, killing
x^k f0. A generating sequence is beta^k mod one lattice L, and two with
one L differ by an automorphism of G, so PW is a function of (a, L). At
M = N, L = Z_2 l + 2^N O0, l primitive, one per point of P^2(Z/2^N)
(ltable()). The characters of G into Z/2^N are c_2(y .), y in Lambda =
{c_2(y l) = 0}, the forms p_k, q_k = c_2(b beta^k) on a basis b1, b2 of
Lambda give G's coordinates, and the kernels of the y off 2 Lambda are
the <h>, h not in 2G, one per point of P^1(Z/2^N). At M < N, Lambda =
<c> + 2^(N-M) P, P = l^perp, c a cyclic summand of P (Smith form), c
over P^1(Z/2^N) on b1, b2 (each Lambda read more than once), c' the
other basis vector, and g_k = (c_2(c beta^k), c_2(c' beta^k) mod 2^M),
whose kernel on (Z/2^N)^3 is L by linearity. A point (s, t) is the
functional (s, 2^(N-M) t) on G, so its kernel is a function of (s, t,
N, M) alone, free of e; the points (1, t), t mod 2^M, and (2 s, 1), s
mod 2^(N-1), reach every <h>, h not in 2G. The character chi with
kernel <h> embeds the cyclic G/<h> in Z/2^N: g_k lies in <h> iff
chi(g_k) = 0, q g_k iff chi(g_k) is a multiple of 2^(N-e), so <h> is
closed iff no active k (q g_k != 0: min of v(F_k), v(2^(N-M) F'_k) at
most N - e - 1) has chi(g_k) a nonzero multiple of 2^(N-e). THE WINDOW:
every a_i is even, so beta^3 lies in 2 O0, beta^(3m) in 2^m O0, and q
beta^k in 2^N O0 from k = 3(N - e): k < 3N reads every active k. At
e = 1, beta a uniformizer, 2 beta^k lies in 2^N O0 exactly from
k = 3N - 3, at every a.

PREDICTIONS (frozen before code; e = 1 first, then e = 2, 3 in a later
run, transplants re-derived above, except where marked added). a =
(2^e, a1, a2), a_i even mod 2^N.
  Q1 (control, exhaustive N = 2, 3, every M in 2..N, every a, a0 = 2
     and 4) every instance's g_k generates G and has l in its kernel;
     the kernels of the instances equal, as a set, the kernels of all
     generating triples of G (both read as masks on (Z/2^N)^3).
  Q2 (control) the engine's closed kernels equal brute_m's closed <h>
     on g_k at every instance, the count equal: e = 1 exhaustive N = 2,
     3, drawn N = 4, 1000 per M, --full N = 5, 500 per M (seed 17321);
     e = 2 exhaustive (8, 8), drawn (16, 8), (16, 16) 1000 each, --full
     N = 5, M = 3..5, 500 each; e = 3 drawn (16, 16) 1000, --full (32,
     16), (32, 32) 500 each (seed 17331, the default's draws first; the
     e = 3 arms added after an audit, since K reads e = 3).
  Q2' (plant, must miss somewhere) the window cut before the last index
     active at that a (3N - 4 at e = 1): off at some instance of every
     drawn arm, not all. Amended at e = 2 after its first run
     (disclosed): cut at 3(N - e) - 1 it was off nowhere at (16, 16),
     since the window bound is loose above e = 1. Q2'' (plant) c'
     unscaled in the activity rule: off at some, not all, below the
     diagonal (at M = N it changes nothing). Q2e (plant, e >= 2) the
     threshold and hit set of e = 1: off at some, not all, on every
     drawn arm. At (32, 32), e = 2, an earlier run (the same seed, its
     draws in another order) found it off at 500 of 500 (its bar failed there); on the diagonal it is
     off at nearly every draw, so it reads little there. A plant off
     nowhere or everywhere reads nothing.
  Q3 (the reader's no) a bare list in (Z/8)^2 found by search with no
     closed <h>, at q = 2 and at q = 4: brute_m and the engine's point
     reader both say none.
  C  (positive controls, known yes) T4 (dualrs.py): a1 = a2 = 0 closed
     at every instance swept; R14: N = M = e + 1 closed at every (a, l),
     read off K.
  W  per arm the least K with q beta^k in 2^N O0 for k >= K, max over a:
     at most 3(N - e); at e = 1, 3N - 3 at every a.
  K  (the kill) per (e, N, M) the count of instances with no closed
     kernel. e = 1: homocyclic N = 2..5 (--full 6), off N = 3, 4 (--full
     5), M = 2..N-1; e = 2: (8, 8), (16, 8), (16, 16), (32, 32), --full
     (32, 8), (32, 16), (64, 64); e = 3: (16, 16), --full (32, 16), (32,
     32). PW at exponent 2^N dies at any nonzero count.
  S  (census) per arm the least closed count in distinct kernels, with
     its (a, l, c).
  T  the group type of O/I at each of dualmult.py's early pairs (j below
     depth) at a q = 4 cubic base with local degree 3: their exponents
     (the bar, all 510 at most 32 and so inside K's reach, added after
     an earlier run printed them).

FINDINGS (entered after the run, from its printed output; --full).
  Q3 q = 2: (2, 0), (0, 2), (2, 2), (2, 4), (4, 6), (3, 0), (6, 2); q =
     4: (7, 3), (0, 0), (3, 1), (5, 0), (5, 4), (4, 5), (2, 5): no closed
     <h> by either reader.
  Q1 at a0 = 2 and 4 alike: 0 not generating, 0 with l off the kernel,
     0 a whose kernels differ: 112 instances at (4, 4), 28 lattices;
     21504 at (8, 4), 168; 1792 at (8, 8), 112.
  Q2 0 off, 0 count off at every arm. e = 1 equals the q = 2 record
     draw for draw: Q2' off at 80, 7168, 1164 exhaustive ((4, 4), (8,
     4), (8, 8)), 111, 325, 655 (N = 4, M = 2..4), 35, 58, 138, 303 (N =
     5, M = 2..5); Q2'' 16032 at (8, 4), 664, 747, 325, 346, 371 below
     the diagonal, 0 on it. e = 2: (8, 8) exhaustive Q2' 1408, Q2e 1764
     of 1792; drawn Q2' 522, 829 (N = 4, M = 3, 4), 123, 196, 373 (N =
     5, M = 3..5), Q2'' 417, 309, 302 below the diagonal, Q2e 921, 394,
     487 below it and 997 of 1000, 497 of 500 on it. e = 3: Q2' 836,
     266, 402 at (16, 16), (32, 16), (32, 32), Q2'' 208 and Q2e 478 at
     (32, 16), Q2e 999 of 1000 and 493 of 500 on the diagonal.
  C  T4 0 fails at every arm; R14's N = M = e + 1, read off K: closed at
     all 112 (q = 2), 1792 (q = 4), 28672 (q = 8).
  W  e = 1: 3N - 3 at every a, N = 2..6. e = 2: max 3, 5, 7, 9 at N =
     3..6 against 3, 6, 9, 12; e = 3: 3, 5 at N = 4, 5 against 3, 6.
  K  0 without a closed <h>. e = 1: M = N at 112, 1792, 28672, 458752,
     7340032 (N = 2..6); 21504 at (8, 4); 688128 at each of (16, 4),
     (16, 8); 22020096 at each of (32, 4), (32, 8), (32, 16). e = 2:
     1792 (8, 8), 688128 (16, 8), 28672 (16, 16), 22020096 each (32, 8),
     (32, 16), 458752 (32, 32), 7340032 (64, 64). e = 3: 28672 (16, 16),
     22020096 (32, 16), 458752 (32, 32). Every case-B type of exponent
     at most 32 at each e: PW, so GC, at q = 2, 4, 8, every cubic
     factor, every G of exponent at most 32, and PW on (Z/64)^2 at q =
     2, 4 (a rule).
  S  least distinct closed kernels. e = 1: 3 at (4, 4), (16, 16); 2 at
     (8, 8), (32, 32), (64, 64) and below the diagonal; 1 at (32, 16), a
     = (2, 2, 2), l = (1, 8, 4), c index 32. e = 2: 3 at (8, 8); 2 at
     (16, 8), (16, 16), (32, 8), (32, 16); 1 at (32, 32), a = (4, 2,
     14), l = (1, 4, 0); 5 at (64, 64). e = 3: 3 at (16, 16), 2 at (32,
     16), (32, 32).
  T  510 early pairs over 18 bases, by the exponent of O/I 4: 130, 8:
     202, 16: 94, 32: 84; every one at most 32. Rank-2 types with both
     factors above 4: (8, 8) 32, (16, 8) 8; the rest are cyclic, rank 3
     or hold a factor at most 4.

RUN RECORD. Default 67/67, 48.3 s, 68.3 MB; --full 119/119, 728.1 s,
before Q2e's four diagonal checks were restored (their values on that
print, 997, 497, 999, 493, pass the bar; not rerun), peak commit
132.4 MB (N = 6 homocyclic 90.9 s at e = 1, 90.4 s at e =
2; the (32, 2^M) arms 63 to 93 s each), both under a memory guard. The
q = 2 script before this one: default 29/29, 25.5 s; --full 48/48,
273.9 s, 110.9 MB, its prints equal to e = 1's here.
"""

import os
import random
import sys
import time
from itertools import combinations
from math import gcd

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

CHECKS = []


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


def omul(x, y, a, M):
    """x y in O0 mod M, beta^3 = -(a2 beta^2 + a1 beta + a0)."""
    c = [0] * 5
    for i in range(3):
        for j in range(3):
            c[i + j] += x[i] * y[j]
    for d in (4, 3):
        t = c[d]
        c[d] = 0
        for i in range(3):
            c[d - 3 + i] -= t * a[i]
    return tuple(v % M for v in c[:3])


def vval(x, M):
    """v_2 elementwise on an int64 array mod M (M's exponent at 0)."""
    x = x % M
    out = np.full(x.shape, M.bit_length() - 1, dtype=np.int64)
    for b in range(M.bit_length() - 2, -1, -1):
        out[(x & ((1 << (b + 1)) - 1)) == (1 << b)] = b
    return out


def vinv(x, M):
    """Inverse of odd x mod M, elementwise (Newton)."""
    y = x % M
    for _ in range(6):
        y = y * (2 - x * y) % M
    return y


def ltable(N):
    """P^2(Z/2^N), first odd coordinate 1, as three int64 arrays."""
    M = 1 << N
    r, e = np.arange(M, dtype=np.int64), np.arange(0, M, 2, dtype=np.int64)
    out = []
    for grid, lead in (((r, r), 0), ((e, r), 1), ((e, e), 2)):
        x, y = [g.ravel() for g in np.meshgrid(*grid, indexing="ij")]
        one = np.ones_like(x)
        out.append((one, x, y) if lead == 0 else
                   (x, one, y) if lead == 1 else (x, y, one))
    return tuple(np.concatenate([o[i] for o in out]) for i in range(3))


def mpoints(N, Mx):
    """The points (s, t) on Lambda of type (2^N, 2^Mx); at Mx = N,
    P^1(Z/2^N)."""
    s = [1] * (1 << Mx) + [2 * x for x in range(1 << (N - 1))]
    t = list(range(1 << Mx)) + [1] * (1 << (N - 1))
    return np.array(s, dtype=np.int64), np.array(t, dtype=np.int64)


def forms(a, l, N):
    """Per l (three arrays), the forms (p_k, q_k) = c_2(b_i beta^k),
    k < 3N, on a basis b1, b2 of Lambda = {y : c_2(y l) = 0 mod 2^N}."""
    M = 1 << N
    pw = [(1, 0, 0)]
    for _ in range(3 * N + 2):
        pw.append(omul(pw[-1], (0, 1, 0), a, M))
    t = [p[2] for p in pw]
    w = [sum(l[x] * t[x + i] for x in range(3)) % M for i in range(3)]
    i0 = np.where(w[0] & 1 == 1, 0, np.where(w[1] & 1 == 1, 1, 2))
    inv = vinv(np.choose(i0, w), M)
    basis = []
    for j in (0, 1):
        # the j-th basis vector: 1 at the j-th index other than i0
        jj = np.where(i0 <= j, j + 1, j)
        wj = np.choose(jj, w)
        basis.append([np.where(i0 == x, -wj * inv % M,
                               np.where(jj == x, 1, 0)) for x in range(3)])
    return [tuple(sum(b[x] * t[x + k] for x in range(3)) % M
                  for b in basis) for k in range(3 * N)]


def coords(fm, N, Mx):
    """G's coordinates (F_k, F'_k) per instance: (p_k, q_k) at Mx = N,
    shape (nl,); else c = al b1 + be b2 over P^1(Z/2^N), c' = b2 where al
    is odd, else b1, shape (nl, nc)."""
    if Mx == N:
        return fm
    M = 1 << N
    al, be = mpoints(N, N)
    odd = (al & 1 == 1)[None, :]
    return [((p[:, None] * al + q[:, None] * be) % M,
             np.where(odd, q[:, None], p[:, None])) for p, q in fm]


def closed(cs, N, Mx, e, kmax=None, unscaled=False, rule_e=None):
    """Per instance and point, whether the point is closed at q = 2^e;
    rule_e reads the threshold and hit set of another e (a plant)."""
    M = 1 << N
    sh = 1 << (N - Mx)
    r = e if rule_e is None else rule_e
    unit = 1 << (N - r)
    s, t = mpoints(N, Mx)
    bad = None
    for F, Fp in cs[:kmax]:
        act = np.minimum(vval(F, M), vval((1 if unscaled else sh) * Fp, M)) \
            <= N - r - 1
        val = (F[..., None] * s + (sh * Fp)[..., None] * t) % M
        hit = (val % unit == 0) & (val != 0) & act[..., None]
        bad = hit if bad is None else bad | hit
    return ~bad


def kernels(N, Mx):
    """Per point, its kernel on Z/2^N x Z/2^Mx, and its kernel id."""
    M, sh = 1 << N, 1 << (N - Mx)
    out, ids = [], {}
    for s, t in zip(*mpoints(N, Mx)):
        K = frozenset((x, y) for x in range(M) for y in range(1 << Mx)
                      if (int(s) * x + sh * int(t) * y) % M == 0)
        out.append(K)
        ids.setdefault(K, len(ids))
    return out, np.array([ids[K] for K in out])


def kcount(cl, kid):
    """Distinct closed kernels per instance."""
    return sum(cl[..., kid == g].any(axis=-1).astype(np.int64)
               for g in range(kid.max() + 1))


def brute_m(g, N, Mx, q):
    """The closed cyclic <h>, h not in 2G, of Z/2^N x Z/2^Mx, by
    membership, at the k with q g_k != 0."""
    M, Q = 1 << N, 1 << Mx
    act = [(x, y) for x, y in g if (q * x) % M or (q * y) % Q]
    out = set()
    for h1 in range(M):
        for h2 in range(Q):
            if not (h1 & 1 or h2 & 1):
                continue
            C = frozenset(((i * h1) % M, (i * h2) % Q) for i in range(M))
            if C in out:
                continue
            if all(((x, y) in C) or (((q * x) % M, (q * y) % Q) not in C)
                   for x, y in act):
                out.add(C)
    return out


def masks(g0, g1, g2, N, Mx):
    """Kernels of x -> sum x_j g_j on (Z/2^N)^3, packed, one row per
    triple; g_j = (F, F') as int64 arrays."""
    M = 1 << N
    r = np.arange(M, dtype=np.int64)
    X = [v.ravel() for v in np.meshgrid(r, r, r, indexing="ij")]
    u = sum(X[j][None, :] * g[0][:, None] for j, g in
            enumerate((g0, g1, g2))) % M
    v = sum(X[j][None, :] * g[1][:, None] for j, g in
            enumerate((g0, g1, g2))) % (1 << Mx)
    return np.packbits((u == 0) & (v == 0), axis=1)


def generating(g0, g1, g2):
    """The mod-2 2x3 matrix (g0 g1 g2) has rank 2."""
    m = [((x[0] & 1), (x[1] & 1)) for x in (g0, g1, g2)]
    out = np.zeros(m[0][0].shape, dtype=bool)
    for i in range(3):
        for j in range(i + 1, 3):
            out |= (m[i][0] * m[j][1] + m[i][1] * m[j][0]) & 1 == 1
    return out


def triple_kernels(N, Mx):
    """The set of kernels of every generating triple of G."""
    els = [(x, y) for x in range(1 << N) for y in range(1 << Mx)]
    E = np.array(els, dtype=np.int64)
    n = len(els)
    out = set()
    for i0 in range(n):
        i1, i2 = [v.ravel() for v in np.meshgrid(np.arange(n), np.arange(n),
                                                 indexing="ij")]
        g0 = (np.full(n * n, E[i0, 0]), np.full(n * n, E[i0, 1]))
        g1, g2 = (E[i1, 0], E[i1, 1]), (E[i2, 0], E[i2, 1])
        keep = generating(g0, g1, g2)
        if not keep.any():
            continue
        sel = lambda g: (g[0][keep], g[1][keep])  # noqa: E731
        out.update(bytes(r) for r in masks(sel(g0), sel(g1), sel(g2), N, Mx))
    return out


def instance_seq(cs, idx, N, Mx):
    """One instance's g_k, k < 3N, as integer pairs."""
    return [(int(F[idx]) % (1 << N), int(Fp[idx]) % (1 << Mx))
            for F, Fp in cs]


def last_active(a, N, e):
    """W: the least K with 2^e beta^k in 2^N O0 for every k >= K."""
    M = 1 << N
    pw = [(1, 0, 0)]
    for _ in range(3 * N + 2):
        pw.append(omul(pw[-1], (0, 1, 0), a, M))
    nz = [k for k, p in enumerate(pw) if any((x << e) % M for x in p)]
    assert not any(any((x << e) % M for x in p) for p in pw[-3:])
    return nz[-1] + 1 if nz else 0


def section_q1():
    print("Q1 the instances' sequences: generating, l in the kernel, every "
          "lattice", flush=True)
    for a0 in (2, 4):
        for N in (2, 3):
            M = 1 << N
            L = ltable(N)
            for Mx in range(2, N + 1):
                t0 = time.time()
                want = triple_kernels(N, Mx)
                ngen = nl_off = ncov = n = 0
                for a1 in range(0, M, 2):
                    for a2 in range(0, M, 2):
                        cs = coords(forms((a0, a1, a2), L, N), N, Mx)
                        shape = cs[0][0].shape
                        lb = [np.broadcast_to(L[x].reshape(
                            (-1,) + (1,) * (len(shape) - 1)), shape)
                            for x in range(3)]
                        n += int(np.prod(shape))
                        ngen += int((~generating(*cs[:3])).sum())
                        u = sum(lb[j] * cs[j][0] for j in range(3)) % M
                        v = sum(lb[j] * cs[j][1] for j in range(3)) % (1 << Mx)
                        nl_off += int(((u != 0) | (v != 0)).sum())
                        flat = [(F.ravel(), Fp.ravel()) for F, Fp in cs[:3]]
                        got = {bytes(r) for r in masks(*flat, N, Mx)}
                        ncov += got != want
                print(f"  a0 {a0}, N {N}, M {Mx}: {n} instances; not "
                      f"generating {ngen}; l off the kernel {nl_off}; a whose "
                      f"kernels differ from the {len(want)} of the generating "
                      f"triples {ncov}; {time.time() - t0:.1f} s", flush=True)
                check(f"Q1 a0 = {a0}, N = {N}, M = {Mx}", ngen == 0 and
                      nl_off == 0 and ncov == 0, f"{ngen}, {nl_off}, {ncov}")


def q2_plants(a, N, e):
    """The plants at a: (name, closed()'s keywords)."""
    out = [("Q2'", {"kmax": last_active(a, N, e) - 1}),
           ("Q2''", {"unscaled": True})]
    if e > 1:
        out.append(("Q2e", {"rule_e": 1}))
    return out


def compare(cs, idx, N, Mx, e, kern, plants):
    """Q2 at one instance: [off, count off] + [plant off per plant]; cs
    sliced to that instance by idx."""
    g = instance_seq(cs, idx, N, Mx)
    br = brute_m(g, N, Mx, 1 << e)
    one = [(np.array([F[idx]]), np.array([Fp[idx]])) for F, Fp in cs]
    got = []
    for kw in ({},) + plants:
        cl = closed(one, N, Mx, e, **kw)[0]
        got.append({kern[i] for i in np.nonzero(cl)[0]})
    return [got[0] != br, len(got[0]) != len(br)] + [x != br for x in got[1:]]


def q2_drawn(rng, e, N, Mx, per):
    """Q2 on per draws of (a, l, c) at one arm."""
    M = 1 << N
    L = ltable(N)
    names = [x for x, _ in q2_plants((1 << e, 0, 0), N, e)]
    kern, _ = kernels(N, Mx)
    tot = np.zeros(2 + len(names), dtype=np.int64)
    for _ in range(per):
        a = (1 << e, 2 * rng.randrange(M // 2), 2 * rng.randrange(M // 2))
        _, plants = zip(*q2_plants(a, N, e))
        i = rng.randrange(len(L[0]))
        la = tuple(np.array([int(x[i])], dtype=np.int64) for x in L)
        cs = coords(forms(a, la, N), N, Mx)
        idx = (0,) if Mx == N else (0, rng.randrange(cs[0][0].shape[1]))
        tot += compare(cs, idx, N, Mx, e, kern, plants)
    pl = "; ".join(f"{x} off {y}" for x, y in zip(names, tot[2:]))
    print(f"  e {e}, N {N}, M {Mx} drawn: {per} instances; off {tot[0]}, "
          f"count off {tot[1]}; {pl}", flush=True)
    check(f"Q2 e = {e}, N = {N}, M = {Mx}", tot[0] == 0 and tot[1] == 0,
          f"{tot[0]}, {tot[1]} of {per}")
    for x, y in zip(names, tot[2:]):
        if x == "Q2''" and Mx == N:
            continue
        check(f"{x} e = {e}, N = {N}, M = {Mx}: off at some, not all",
              0 < y < per, f"{y} of {per}")


def section_q2(full):
    print("Q2 the engine against brute_m (plants: Q2' window, Q2'' "
          "unscaled, Q2e the e = 1 rule)", flush=True)
    for e, arms in ((1, [(2, 2), (3, 2), (3, 3)]), (2, [(3, 3)])):
        for N, Mx in arms:
            M = 1 << N
            L = ltable(N)
            names = [x for x, _ in q2_plants((1 << e, 0, 0), N, e)]
            kern, _ = kernels(N, Mx)
            tot = np.zeros(2 + len(names), dtype=np.int64)
            n = 0
            for a1 in range(0, M, 2):
                for a2 in range(0, M, 2):
                    _, plants = zip(*q2_plants((1 << e, a1, a2), N, e))
                    cs = coords(forms((1 << e, a1, a2), L, N), N, Mx)
                    for idx in np.ndindex(cs[0][0].shape):
                        tot += compare(cs, idx, N, Mx, e, kern, plants)
                        n += 1
            pl = "; ".join(f"{x} off {y}" for x, y in zip(names, tot[2:]))
            print(f"  e {e}, N {N}, M {Mx} exhaustive: {n} instances; off "
                  f"{tot[0]}, count off {tot[1]}; {pl}", flush=True)
            check(f"Q2 e = {e}, N = {N}, M = {Mx}", tot[0] == 0 and
                  tot[1] == 0, f"{tot[0]}, {tot[1]} of {n}")
    rng = random.Random(17321)
    for N, per in ((4, 1000),) + (((5, 500),) if full else ()):
        for Mx in range(2, N + 1):
            q2_drawn(rng, 1, N, Mx, per)
    rng = random.Random(17331)
    drawn = [(2, 4, 3, 1000), (2, 4, 4, 1000), (3, 4, 4, 1000)]
    if full:
        drawn += [(2, 5, Mx, 500) for Mx in range(3, 6)] + \
            [(3, 5, 4, 500), (3, 5, 5, 500)]
    for e, N, Mx, per in drawn:
        q2_drawn(rng, e, N, Mx, per)


def section_q3():
    print("Q3 the readers' no: a bare list in (Z/8)^2", flush=True)
    for e, seed in ((1, 17322), (2, 17332)):
        rng = random.Random(seed)
        for _ in range(200000):
            lst = [(rng.randrange(8), rng.randrange(8))
                   for _ in range(rng.randrange(3, 9))]
            if not brute_m(lst, 3, 3, 1 << e):
                break
        one = [(np.array([x]), np.array([y])) for x, y in lst]
        n = int(closed(one, 3, 3, e)[0].sum())
        check(f"Q3 q = {1 << e}: brute_m and the point reader both find no "
              f"closed <h>", not brute_m(lst, 3, 3, 1 << e) and n == 0,
              f"{lst}, engine closed {n}")


def section_k(full):
    print("K, S, C, W the sweep: instances with no closed kernel, the least "
          "closed count, the window", flush=True)
    arms = [(1, N, N) for N in (2, 3, 4, 5)] + \
        [(1, N, Mx) for N in (3, 4) for Mx in range(2, N)] + \
        [(2, 3, 3), (2, 4, 3), (2, 4, 4), (2, 5, 5), (3, 4, 4)]
    if full:
        arms += [(1, 6, 6)] + [(1, 5, Mx) for Mx in range(2, 5)] + \
            [(2, 5, 3), (2, 5, 4), (2, 6, 6), (3, 5, 4), (3, 5, 5)]
    for e, N, Mx in arms:
        t0 = time.time()
        M = 1 << N
        L = ltable(N)
        _, kid = kernels(N, Mx)
        none = t4 = n = 0
        wmin, wmax = 3 * N, 0
        least = None
        for a1 in range(0, M, 2):
            for a2 in range(0, M, 2):
                a = (1 << e, a1, a2)
                w = last_active(a, N, e)
                wmin, wmax = min(wmin, w), max(wmax, w)
                cnt = kcount(closed(coords(forms(a, L, N), N, Mx), N, Mx, e),
                             kid)
                n += cnt.size
                none += int((cnt == 0).sum())
                if a1 == 0 and a2 == 0:
                    t4 += int((cnt == 0).sum())
                idx = np.unravel_index(np.argmin(cnt), cnt.shape)
                if least is None or cnt[idx] < least[0]:
                    least = (int(cnt[idx]), a,
                             tuple(int(x[idx[0]]) for x in L),
                             None if Mx == N else int(idx[1]))
        print(f"  e {e}, N {N}, M {Mx}: {n} instances; K no closed <h> at "
              f"{none}; C T4 fails at {t4}; W window {wmin}..{wmax} over a "
              f"(bound {3 * (N - e)}); S least {least[0]} distinct closed "
              f"kernels at a = {least[1]}, l = {least[2]}, c index "
              f"{least[3]}; {time.time() - t0:.1f} s", flush=True)
        check(f"K PW at e = {e}, N = {N}, M = {Mx}", none == 0,
              f"{none} of {n}")
        check(f"C T4 at e = {e}, N = {N}, M = {Mx}", t4 == 0, f"{t4}")
        check(f"W window at e = {e}, N = {N}, M = {Mx}",
              wmax <= 3 * (N - e) and (e > 1 or wmin == wmax == 3 * N - 3),
              f"{wmin}..{wmax}")


def smith(rows):
    """The invariant factors > 1 of Z^3 / rows, largest first, from the
    determinantal divisors."""
    n = len(rows)
    dets = [1]
    for r in range(1, n + 1):
        g = 0
        for I in combinations(range(n), r):
            for J in combinations(range(n), r):
                m = [[rows[i][j] for j in J] for i in I]
                g = gcd(g, int(round(abs(np.linalg.det(np.array(
                    m, dtype=float))))))
        dets.append(g)
    inv = [dets[i] // dets[i - 1] for i in range(1, n + 1)]
    return tuple(sorted((x for x in inv if x > 1), reverse=True))


def section_t():
    print("T the early pairs at q = 4 cubic bases of local degree 3: the "
          "type of O/I", flush=True)
    import dualmult
    t0 = time.time()
    bases = [a for a in dualmult.cubic_bases()
             if abs(a[0]) == 4 and a[1] % 2 == 0 and a[2] % 2 == 0]
    rows, _ = dualmult.run(bases, lambda q: 256)
    early = [r for r in rows if r["dep"] is not None and r["j"] < r["dep"]]
    types, expo = {}, {}
    for r in early:
        ty = smith(r["L"].rows())
        types[ty] = types.get(ty, 0) + 1
        expo[ty[0]] = expo.get(ty[0], 0) + 1
    print(f"  {len(bases)} bases, {len(early)} early pairs; by type "
          f"{dict(sorted(types.items()))}; by exponent "
          f"{dict(sorted(expo.items()))}; {time.time() - t0:.1f} s",
          flush=True)
    check("T the early pairs at local degree 3 are dualmult.py's 510, each "
          "of exponent at most 32", len(early) == 510 and max(expo) <= 32,
          f"{len(early)}, max {max(expo)}")


def main():
    full = "--full" in sys.argv
    t0 = time.time()
    section_q3()
    section_q1()
    section_q2(full)
    section_k(full)
    section_t()
    print(f"{sum(CHECKS)}/{len(CHECKS)} checks, {time.time() - t0:.1f} s")
    sys.exit(0 if all(CHECKS) else 1)


if __name__ == "__main__":
    main()