import numpy as np
rng=np.random.default_rng(31)
def step(a,d):
    L=len(a); n=np.zeros(L,dtype=np.uint8)
    for k in range(L):
        left=a[k-1] if k>0 else d
        right=a[k+1] if k+1<L else 0
        n[k]=right^(a[k]|left)
    return n
def feats(a):
    N=int(a.sum()); R=int(((a[1:]==1)&(a[:-1]==0)).sum()+(a[0]==1))
    P11=int((a[1:]&a[:-1]).sum()); C=int((a[:-2]&(1-a[1:-1])&a[2:]).sum())
    pos=np.arange(1,len(a)+1); ev=(pos%2==0); st=np.concatenate([[a[0]==1],(a[1:]==1)&(a[:-1]==0)])
    Ne=int(a[ev].sum()); Re=int(st[ev].sum()); P11e=int((a[1:]&a[:-1])[ev[:-1]].sum())
    return [N%2,R%2,P11%2,C%2,Ne%2,Re%2,P11e%2,int(a[0]),int(a[1]),int(a[2])]
names=[]
for tt in ("t","t+1","t+2"):
    names+= [f"N({tt})",f"R({tt})",f"P11({tt})",f"C101({tt})",f"Neven({tt})",f"Reven({tt})",f"P11even({tt})",f"a1({tt})",f"a2({tt})",f"a3({tt})"]
names.append("1")
rows=[]
for _ in range(6000):
    L=int(rng.integers(6,30)); a=np.zeros(L+8,dtype=np.uint8); a[:L]=rng.integers(0,2,size=L)
    f=[]; x=a
    for d in (1,0,1):                     # lit, dark, lit
        f+=feats(x); x=step(x,d)
    rows.append(f+[1])
M=np.array(rows,dtype=np.uint8)%2
# nullspace over GF(2): find all c with M c = 0
A=M.copy(); n=A.shape[1]; piv=[]; r=0
for c in range(n):
    p=next((i for i in range(r,len(A)) if A[i,c]),None)
    if p is None: continue
    A[[r,p]]=A[[p,r]]
    for i in range(len(A)):
        if i!=r and A[i,c]: A[i]^=A[r]
    piv.append(c); r+=1
free=[c for c in range(n) if c not in piv]
basis=[]
for f in free:
    v=np.zeros(n,dtype=np.uint8); v[f]=1
    for i,c in enumerate(piv):
        if A[i,f]: v[c]=1
    basis.append(v)
print(f"exact mod-2 relations holding for ALL 6000 random crowds: {len(basis)}")
for v in basis: print("   ", " + ".join(names[i] for i in range(n) if v[i]), "= 0")