Counting permutations with roots

My publish from yesterday on permutation roots ends with a Mathematica code for locating the likelihood {that a} permutation of n components has a okayth root. That is carried out by discovering the coefficient of xⁿ within the producing operate

I wished to say extra about this, and take a look at implementing the identical code in SymPy. I used to be curious how nicely SymPy would do as a result of I’ve observed that LLMs typically generate SymPy code because it’s an open supply CAS.

Wilf [1] describes the infinite product above because the exponential producing operate (egf) of f(n, okay), the variety of permutations of n objects which have a okayth root. Since egfs have a n! time period within the denominator, that is additionally the odd producing operate (ogf) of the likelihood {that a} randomly chosen permutation on n objects has a okayth root.

My first try at utilizing Mathematica to probe the producing operate was

expq[x_, q_] := MittagLefflerE[q, x^q]	 
p[n_, k_] :=  SeriesCoefficient[	 
    Product[expq[x^m/m, GCD[m, k]], {m, 1, Infinity}], {x, 0, n}]

This hung endlessly after I tried to apply it to a small instance. I noticed, however apparently Mathematica didn’t, that Infinity could possibly be changed by n since phrases increased than n don’t contribute to the coefficient of xⁿ. With that change, the code ran shortly.

This morning I attempted changing the Mathematica code to Sympy; Claude did this in a single shot. I additionally reproduced the desk of f(n, okay) values on web page 150 of [1] to check the code. Since Wilf tabulated f(n, okay), not f(n, okay)/n!, I multiplied the outcomes by n!.

Right here is the output:

okay = 2 [1, 1, 3, 12, 60, 270, 1890, 14280, 128520, 1096200]
okay = 3 [1, 2, 4, 16, 80, 400, 2800, 22400, 181440, 1814400]
okay = 4 [1, 1, 3, 12, 60, 270, 1890, 13020, 117180, 1039500]
okay = 5 [1, 2, 6, 24, 96, 576, 4032, 32256, 290304, 2612736]
okay = 6 [1, 1, 1, 4, 40, 190, 1330, 8680, 52920, 340200]
okay = 7 [1, 2, 6, 24, 120, 720, 4320, 34560, 311040, 3110400]

and right here is the SymPy code. I edited the principle however the remaining is verbatim from Claude.

from sympy import symbols, gcd, factorial, Rational, S

x = symbols('x')

def expq_coeffs(m, q, n):
    """
    Truncated (diploma <= n) sequence coefficients of
        expq(x**m/m, q) = MittagLefflerE(q, (x**m/m)**q)
    Since q is a optimistic integer:
        E_q(y^q) = sum_j y^(q*j) / (q*j)!
    with y = x**m/m, so the time period of diploma m*q*j has coefficient
        1 / ( m**(q*j) * (q*j)! ).
    Returns an inventory c[0..n] of coefficients.
    """
    c = [S.Zero] * (n + 1)
    j = 0
    whereas m * q * j <= n:
        deg = m * q * j
        c[deg] += Rational(1, m**(q * j) * factorial(q * j))
        j += 1
    return c

def poly_mult_trunc(a, b, n):
    """Multiply two sequence (lists of coeffs, index = diploma) truncated to diploma n."""
    c = [S.Zero] * (n + 1)
    for i, ai in enumerate(a):
        if ai == 0:
            proceed
        max_j = n - i
        for j2 in vary(max_j + 1):
            bj = b[j2]
            if bj != 0:
                c[i + j2] += ai * bj
    return c

def p(n, okay):
    """
    SymPy equal of:
        expq[x_, q_] := MittagLefflerE[q, x^q]
        p[n_, k_] := SeriesCoefficient[
            Product[expq[x^m/m, GCD[m, k]], {m, 1, n}], {x, 0, n}]
    """
    outcome = [S.Zero] * (n + 1)
    outcome[0] = S.One
    for m in vary(1, n + 1):
        q = gcd(m, okay)
        issue = expq_coeffs(m, q, n)
        outcome = poly_mult_trunc(outcome, issue, n)
    return outcome[n]

# instance
if __name__ == "__main__":
    for okay in vary(2, 8):
        print("okay =", okay, [factorial(n)*p(n, k) for n in range(1,11)])

[1] Herbert Wilf. Generatingfunctionology. Accessible on-line right here.