Counting permutations with roots

(832) 422-8646

Contact

My post from yesterday on permutation roots ends with a Mathematica code for finding the probability that a permutation of n elements has a kth root. This is done by finding the coefficient of xn in the generating function

I wanted to say more about this, and look at implementing the same code in SymPy. I was curious how well SymPy would do because I’ve noticed that LLMs often generate SymPy code since it’s an open source CAS.

Wilf [1] describes the infinite product above as the exponential generating function (egf) of f(n, k), the number of permutations of n objects that have a kth root. Since egfs have a n! term in the denominator, this is also the ordinary generating function (ogf) of the probability that a randomly chosen permutation on n objects has a kth root.

My first attempt at using Mathematica to probe the generating function was

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

This hung forever when I tried to use it on a small example. I realized, but apparently Mathematica did not, that Infinity could be replaced by n since terms higher than n do not contribute to the coefficient of xn. With that change, the code ran quickly.

This morning I tried converting the Mathematica code to Sympy; Claude did this in one shot. I also reproduced the table of f(n, k) values on page 150 of [1] to test the code. Since Wilf tabulated f(n, k), not f(n, k)/n!, I multiplied the results by n!.

Here is the output:

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

and here is the SymPy code. I edited the main but the rest is verbatim from Claude.

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

x = symbols('x')

def expq_coeffs(m, q, n):<br>"""<br>Truncated (degree<br>[1] Herbert Wilf. Generatingfunctionology. Available online here.

Leave a Reply<br>Your email address will not be published. Required fields are marked *<br>Comment *<br>Name *

Email *

Website

Search for:

John D. Cook, PhD

My colleagues and I have decades of consulting experience helping companies solve complex problems involving data privacy, applied math, and statistics.

Let’s talk. We look forward to exploring the opportunity to help your company too.