def fdeg(p):

    R.<q> = 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.<q> = 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.<q> = PolynomialRing(QQ)

    R.<q> = 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

​

​

parent_initializer = SkewPartitions

element_repr = lambda P: str(list(P)).replace(" ","")

levels = [2, 3, 4, 5]

for level in levels:

    for elt in parent_initializer(level):

        print('%s  =>  %s' % (element_repr(elt), statistic(elt)))