How many elements of order $k$ are in $S_n$?

Consider all cycle types of permutations in $S_{n}$ which have order $k$. Recalling that $\sigma, \tau \in S_{n}$ are conjugate if and only if they have the same cycle type, then it should just then be a matter of finding the number of conjugates of each cycle type with order $k$. This is given by

$$\frac{n!}{\prod_{i=1}^{s}(k_{i}! m_{i}^{k_{i}})}$$

where for each $i \in \{0, 1, \dots, s\}$, $m_i$ is a distinct integer appearing in the cycle type of your permutation, and $k_i$ gives the number of cycles of length $m_{i}$.

This one can be done by inclusion-exclusion. Note that the combinatorial class $\mathcal{Q}$ of permutations whose order divides $k$ is given by $$\def\textsc#1{\dosc#1\csod} \def\dosc#1#2\csod{{\rm #1{\small #2}}}\mathcal{Q} = \textsc{SET}\left(\sum_{d|k} \textsc{CYC}_{=d}(\mathcal{Z})\right).$$ Translation to exponential generating functions we obtain the EGF of permutations whose order divides $k$, which is $$Q_k(z) = \exp\left(\sum_{d|k} \frac{z^d}{d}\right).$$ Now we can use this generating function to count permutations of order exactly $k$ by inclusion-exclusion. Consider the Hasse diagram of the divisor poset of $k.$ By inclusion-exclusion we label the top node, which is $k$, with a weight of one representing permutations of order exactly $k$ and descend along the chains, labelling every node $d$ with the appropriate weight to make the weights of the poset spanned by $d$ and $k$ sum to zero.

Let $f_1(d)$ be the weight function. The above process gives the following two equations: $$f_1(k) = 1 \quad\text{and for}\quad d|k,\; d<k,\quad \sum_{m|k/d} f_1(dm) = 0.$$ Now introduce $f_2(d) = f_1(k/d).$ The defining equations then become $$f_2(1) = 1 \quad\text{and for}\quad d|k,\; d>1,\quad \sum_{m|d} f_2(m) = 0.$$ But this is the defining equation of the Mobius function $\mu(d)$ from basic number theory. Therefore the weight of the inclusion-exclusion terms is given by $\mu(k/d)$ and we finally have the EGF $$Q(z) = \sum_{d|k} \mu(k/d) \times Q_d(z) = \sum_{d|k} \mu(k/d) \exp\left(\sum_{m|d} \frac{z^m}{m}\right).$$ The desired count is then given by $$n! [z^n] Q(z).$$

This formula produces e.g. for $k=6$ the EGF $$Q(z) = {{\rm e}^{z}}-{{\rm e}^{z+1/2\,{z}^{2}}}-{{\rm e}^{z+1/3\,{z}^{3}}}+{{\rm e}^{z+1/2\,{ z}^{2}+1/3\,{z}^{3}+1/6\,{z}^{6}}}$$ with the sequence of values starting at $n=5$ $$20, 240, 1470, 10640, 83160, 584640, 4496030, 42658440, 371762820, 3594871280,\ldots$$ which is OEIS A061121.

For $k=8$ we get the EGF $$Q(z) = -{{\rm e}^{z+1/2\,{z}^{2}+1/4\,{z}^{4}}}+{{\rm e}^{z+1/2\,{z}^{2}+1/4\,{z}^{4}+1/8\,{z }^{8}}}$$ with the sequence of values starting at $n=8$ $$5040, 45360, 453600, 3326400, 39916800, 363242880, 3874590720, 34767532800,\ldots$$ which is OEIS A061122.

Finally for $k=12$ we get the EGF $$Q(z) = {{\rm e}^{z+1/2\,{z}^{2}}}-{{\rm e}^{z+1/2\,{z}^{2}+1/4\,{z}^{4}}}-{{\rm e}^{z+1/2\,{z }^{2}+1/3\,{z}^{3}+1/6\,{z}^{6}}}+{{\rm e}^{z+1/2\,{z}^{2}+1/3\,{z}^{3}+1/4\,{z}^{4}+1 /6\,{z}^{6}+1/12\,{z}^{12}}}$$ with the sequence of values starting at $n=7$ $$420, 3360, 30240, 403200, 4019400, 80166240, 965284320, 12173441280, 162850287600,\ldots$$ which is OEIS A061125.

Remark as of Mon Jul 28 2014. The argument above is correct but it can be simplified by recognizing the Moebius inversion formula that appears.

Addendum as of Sat Apr 22 2017. As rightly pointed out by @Nesa the EGF yields the correct data but is impracticable for $n$ large. We may instead use the recursive formula by Lovasz for the cycle index of the symmetric group, which is

$$Z(S_n) = \frac{1}{n} \sum_{l=1}^n a_l Z(S_{n-l}) \quad\text{where}\quad Z(S_0) = 1.$$

Restricting to cycles of length a divisor of $k$ we obtain

$$Z_k(S_n) = \frac{1}{n} \sum_{l|k \wedge 1\le l\le n} a_l Z_k(S_{n-l}) \quad\text{where}\quad Z_k(S_0) = 1.$$

We then retain only those conjugacy classes (which are represented by partitions) with LCM of the cycle lengths equal to $k$. Further optimization is possible here (we tried to keep it simple) which the reader is invited to investigate. Some examples for the optimized routine compared to Taylor series applied to the EGF and implemented in Maple are

• permutations on $50$ elements of order exactly $24$ in $0.8s$ (as opposed to $31s$, $831$ Mb memory allocated):

$$62691738121269819754220487957220800432035041174469549752320000$$

• permutations on $50$ elements of order exactly $30$ in $0.6s$ (as opposed to $35s$, $846$ Mb memory allocated):

$$105557285434280619538495502609995618378421158878971545429804800$$

• permutations on $50$ elements of order exactly $35$ in $0.01s$ (as opposed to $13s$, $595$ Mb memory allocated):

$$14311697149323256258802803652030609474009628711255472000512000.$$

The Maple code that was used to compare the EGF to the cycle index goes as follows.

with(numtheory):

pet_cycleind_symm_divk :=
proc(n, k)
option remember;
local lseq;

    if n=0 then return 1; fi;

    lseq := select(l -> l <= n, divisors(k));

    expand(1/n*
           add(a[l]*
               pet_cycleind_symm_divk(n-l, k), l in lseq));
end;

QAUX := m -> exp(add(z^d/d, d in divisors(m)));

Q := k -> add(mobius(k/d)*QAUX(d), d in divisors(k));

A := (n, k) -> n!*coeftayl(Q(k), z=0, n);

AOPT :=
proc(n, k)
option remember;
local idx, term, res, ord;

    if n = 1 then
        idx := [a[1]]
    else
        idx := pet_cycleind_symm_divk(n, k);
    fi;

    res := 0;
    for term in idx do
        ord :=
        lcm(seq(op(1, var), var in indets(term)));

        if ord = k then
            res := res + lcoeff(term);
        fi;
    od;

    n!*res;
end;