files / problems / covering / compute / q10 / prescribe_group.py
#!/usr/bin/env python3
"""Prescribe a 2-group and look for an invariant radius-2 set of a given size.
q10's odd-prime sweep does not touch p = 2, and aut_group.c says that is
exactly where the known optima live: the certified 50-set has |Aut| = 128 and
the Kaikkonen-Rosendahl 51-set has |Aut| = 512, both 2-groups. So prescribe
those groups, and their subgroups, instead.
An H-invariant set is a union of H-orbits, so its size is a sum of orbit sizes
-- for a 2-group every orbit size is a power of two, and an odd target forces
at least one H-fixed vector. This enumerates the unions of the right total
size and runs the flat syndrome sweep on each.
Usage: prescribe_group.py <listing> <n> [--subgroups DEPTH] [--cap N]
"""
import sys
from itertools import combinations
def load(path):
basis = S = None
gens = []
for line in open(path):
t = line.split()
if not t:
continue
if t[0] == "basis":
basis = [int(x) for x in t[1:]]
elif t[0] == "set":
S = [int(x) for x in t[1:]]
elif t[0] == "g":
gens.append([int(x) for x in t[1:]])
r = len(basis)
def to_cols(img):
piv = []
for j in range(r):
cur_s, cur_d = basis[j], img[j]
for ps, pd in piv:
if (cur_s >> (ps.bit_length() - 1)) & 1:
cur_s ^= ps
cur_d ^= pd
if cur_s:
piv.append((cur_s, cur_d))
cols = []
for i in range(r):
v, out = 1 << i, 0
for ps, pd in piv:
if (v >> (ps.bit_length() - 1)) & 1:
v ^= ps
out ^= pd
assert v == 0
cols.append(out)
return tuple(cols)
return r, S, [to_cols(g) for g in gens]
def apply(cols, v, r):
out = 0
for i in range(r):
if (v >> i) & 1:
out ^= cols[i]
return out
def close(gens, r):
ident = tuple(1 << i for i in range(r))
G, frontier = {ident}, [ident]
while frontier:
a = frontier.pop()
for b in gens:
c = tuple(apply(a, b[i], r) for i in range(r))
if c not in G:
G.add(c)
frontier.append(c)
return G
def orbits(G, r):
N = 1 << r
seen = [False] * N
out = []
for v in range(1, N):
if seen[v]:
continue
orb = set()
stack = [v]
seen[v] = True
while stack:
u = stack.pop()
orb.add(u)
for g in G:
t = apply(g, u, r)
if not seen[t]:
seen[t] = True
stack.append(t)
out.append(sorted(orb))
return out
def covers(cols, r):
N = 1 << r
seen = bytearray(N)
seen[0] = 1
for a in cols:
seen[a] = 1
for i, a in enumerate(cols):
for b in cols[i + 1:]:
seen[a ^ b] = 1
return N - sum(seen)
def search(orbs, n, r, cap):
"""every union of orbits of total size n, flat-tested"""
orbs = sorted(orbs, key=len)
sizes = [len(o) for o in orbs]
suffix = [0] * (len(orbs) + 1)
for i in range(len(orbs) - 1, -1, -1):
suffix[i] = suffix[i + 1] + sizes[i]
best = [1 << 30, None]
tried = [0]
def rec(i, total, chosen):
if tried[0] > cap:
return
if total == n:
tried[0] += 1
cols = [x for j in chosen for x in orbs[j]]
h = covers(cols, r)
if h < best[0]:
best[0], best[1] = h, list(chosen)
return
if total > n or i == len(orbs) or total + suffix[i] < n:
return
rec(i + 1, total + sizes[i], chosen + [i])
rec(i + 1, total, chosen)
rec(0, 0, [])
return tried[0], best
def index2_subgroups(G, r):
"""kernels of the homomorphisms G -> C_2, i.e. the index-2 subgroups"""
Gl = sorted(G)
idx = {g: i for i, g in enumerate(Gl)}
# G/[G,G]G^2 as an F_2 space: build the quotient by brute force on G
subs = []
for k in range(len(Gl)):
# candidate: the set of g with a chosen linear functional 0 -- instead
# just test every subset of index 2 generated by dropping one generator
pass
# practical route: every maximal subgroup of a 2-group has index 2 and
# contains the Frattini subgroup, so enumerate subgroups generated by all
# but one element of a small generating set
return subs
def main():
path, n = sys.argv[1], int(sys.argv[2])
cap = 200000
depth = 0
if "--cap" in sys.argv:
cap = int(sys.argv[sys.argv.index("--cap") + 1])
if "--subgroups" in sys.argv:
depth = int(sys.argv[sys.argv.index("--subgroups") + 1])
r, S, gens = load(path)
G = close(gens, r)
orbs = orbits(G, r)
from collections import Counter
print(f"{path}: |G|={len(G)} orbits on F_2^{r}\\{{0}}: "
f"{sorted(Counter(len(o) for o in orbs).items())}")
tried, best = search(orbs, n, r, cap)
print(f" unions of total size {n}: {tried} tested, best holes = "
f"{best[0] if best[1] is not None else 'none exist'}")
if depth:
# walk down through subgroups generated by proper subsets of a small
# generating set, so the orbits refine and more sizes become reachable
base = sorted(G)
seen = set()
for k in range(1, depth + 1):
for sub in combinations(base[1:min(len(base), 14)], k):
H = close(list(sub), r)
key = frozenset(H)
if key in seen or len(H) == len(G):
continue
seen.add(key)
ho = orbits(H, r)
t, b = search(ho, n, r, cap)
if b[1] is not None:
print(f" subgroup |H|={len(H)} orbits={len(ho)}: "
f"{t} unions tested, best holes = {b[0]}")
if b[0] == 0:
cols = [x for j in b[1] for x in ho[j]]
print(" COVERING FOUND:", " ".join(map(str, sorted(cols))))
return
if __name__ == "__main__":
main()