dualpw.py
Checks: DUAL: The chain at depth theorem; DUAL: The witness at q = 2^e rule.
Cited on: DUAL.
To run it: save dualpw.py, with dualdeg.py, dualmult.py, dualn.py beside it; needs numpy.
"""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()