def restrict(G, domain):
    H_gens = [[x for x in g.cycle_tuples() if x[0] in domain] for g in G.gens()]
    H = PermutationGroup(H_gens, domain=domain)
    # this normalizes the domain to 1..n
    K = PermutationGroup(gap_group=H._gap_(), domain=range(1, len(domain)+1)) 
    return K

def find_factor(G):
    S = SymmetricGroup(G.degree())
    for a, b in sorted(SetPartitions(G.orbits(), 2), key=lambda x: len(x[0])):
        A = restrict(G, sum(a, tuple()))
        B = restrict(G, sum(b, tuple()))
        if libgap.IsConjugate(S, A.direct_product(B, False), G):
            return A, B

def factor(G):
    F = []
    B = G
    while True:
        f = find_factor(B)
        if f is None:
            F.append(B)
            return [(H, 1) for H in F]
        A, B = f
        F.append(A)

def statistic(G):
    H = G.relabel(range(1, len(G)+1), inplace=False).automorphism_group()
    return sum(f[1] for f in factor(H))