Counting permutations with roots

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] p[n_, k_] := SeriesCoefficient[ Product[expq[x^m/m, GCD[m, k]], {m, 1, Infinity}], {x, 0, n}] This hung…