def statistic(pi):
    L = pi.cycle_type()
    k = pi.number_of_fixed_points()
    if pi(1) == 1:
        k -= 1
    m = len(L)    
    n = len(pi)
    return factorial(n+m-2*(k+1))/factorial(n-k)*prod(L)

# alternative code for checking
@cached_function
def number_of_minimal_star_factorizations_dict(n):
    S = Permutations(n)
    g = [S(Permutation((1, i))) for i in range(2, n+1)]

    result = dict()
    result[S(Permutation([]))] = 1
    l = 1 # the length of the minimal factorization
    while len(result) < factorial(n):
        result_l = dict()
        for w in cartesian_product([g]*l):
            pi = prod(w)
            if pi not in result:
                result_l[pi] = result_l.get(pi, 0) + 1
        for pi, v in result_l.items():
            result[pi] = v
        l += 1
    return result

def statistic(pi):
    d = number_of_minimal_star_factorizations_dict(len(pi))
    return d[pi]