def fdeg(p):
    R. = QQ[]
    f = prod([ 1-q*q^j for j in range(p.size()) ])
    b = sum([ j*k for j, k in enumerate(p) ])
    return (q^b)*(f/p.hook_polynomial(q,q))

def fake_degree(a):
    R. = QQ[]
    if a == 0:
        return R(0)
    else:
        z = SymmetricFunctions(QQ).schur()(a)
        return R(sum([z.coefficient(p)*fdeg(p) for p in z.support()]))

def standard_CSP(p, r):
    P = parent(p)
    q = P.gen()
    PP. = PolynomialRing(QQ)
    R. = PP.quo(q^r-1)
    return R(p).lift()

def orbit_lengths(p, r):
    p = standard_CSP(p, r)
    p = {k: p[k] for k in range(r)}
    divs = sorted(divisors(r))
    O = dict()
    for d in divs:
        count = p[r-d]
        if count:
            O[r//d] = count
        if count:
            for k in range(0, r, d):
                p[k] -= count
    assert all(c == 0 for c in p.values()), "%s" % p
    return O

def has_cyclic_permutation_representation(p, r):
    if p == 0:
        return True
    try:
        lo = orbit_lengths(p, r)
    except AssertionError:
        return False
    return all(v in ZZ and ZZ(v) >= 0 for v in lo.values())

def statistic(la):
    s = SymmetricFunctions(QQ).schur()
    f = fake_degree(s(la))
    for r in reversed(divisors(la.size())):
        if has_cyclic_permutation_representation(f, r):
            return r

# alternative implementation using Theorem 2.7
def has_cyclic_permutation_representation_slow(p, r):
    for k in divisors(r):
        if sum(moebius(k//j) * p.subs(q=QQbar.zeta(r)^j) for j in divisors(k)) < 0:
            return False
    return True

def statistic_slow(la):
    s = SymmetricFunctions(QQ).schur()
    f = fake_degree(s(la))
    for r in reversed(divisors(la.size())):
        if has_cyclic_permutation_representation_slow(f, r):
            return r