Probability: Permutations 
Consider the experiment of picking a random permutation $\pi$ on $\{1,2,...,n\}$, and define the associated random variable $f(\pi)$ as the number of fixed points of $\pi$, i.e, the number of $i$ such that $f(i)=i$.

I know that a permutation of $X=\{1,2,\ldots ,n\}$ is a one-to-one function $\pi: X \rightarrow X$. Any two functions, $\pi_1, \pi_2$ can be composed and the resulting function is also a permutation. How can I find $E(F)$ and what do we know about $E(f(\pi^2))$ and $ E(f( \pi^k))$? 
 A: These can be done using generating functions. First, consider $E[F].$ The exponential generating function of the set of permutations by sets of cycles where fixed points are marked is
$$ G(z, u) = \exp\left(uz - z + \log \frac{1}{1-z}\right) =
\frac{1}{1-z} \exp(uz-z).$$
Now to get $E[F]$ compute
$$ \left.\frac{d}{du} G(z, u)\right|_{u=1} =
\left. \frac{u^0 z}{1-z} \exp(uz-z) \right|_{u=1} = \frac{z}{1-z}.$$
This means that $E[F] = 1,$ there is one fixed point on average.
Let $E[F_2]$ denote the expectation of the number of fixed points in $\sigma^2$, where $\sigma$ is a random permutation. Now we need to mark both fixed points and two-cycles, since the latter turn into two fixed points under squaring. We get
$$ G(z, u) = \exp\left(uz - z + u^2 \frac{z^2}{2} -  \frac{z^2}{2}
+ \log \frac{1}{1-z}\right) =
\frac{1}{1-z} \exp\left(uz - z + u^2 \frac{z^2}{2} -  \frac{z^2}{2}\right).$$
Continuing as before, we find
$$ \left.\frac{d}{du} G(z, u)\right|_{u=1} =
\left. \frac{u^0 z+ u^1 z^2}{1-z} 
\exp\left(uz - z + u^2 \frac{z^2}{2} -  \frac{z^2}{2}\right) \right|_{u=1} =
\frac{z+z^2}{1-z}.$$
The conclusion is that $E[F_2] = 2$ for $n\ge 2$ and there are two fixed points on average.
The pattern should now be readily apparent. For every divisor $d$ of $k$ a cycle of length $d$ splits into $d$ fixed points when raised to the power $k.$ Hence we need to mark these cycles with $u^d.$ To illustrate this consider $E[F_6].$
We get
$$ G(z, u) = \exp\left(uz - z + u^2 \frac{z^2}{2} -  \frac{z^2}{2} +
u^3 \frac{z^3}{3} -  \frac{z^3}{3} + u^6 \frac{z^6}{6} -  \frac{z^6}{6} 
+ \log \frac{1}{1-z}\right) \\ =
\frac{1}{1-z} \exp\left(uz - z + u^2 \frac{z^2}{2} -  \frac{z^2}{2} +
u^3 \frac{z^3}{3} -  \frac{z^3}{3} + u^6 \frac{z^6}{6} -  \frac{z^6}{6}\right).$$
Once more continuing as before, we find
$$ \left.\frac{d}{du} G(z, u)\right|_{u=1} \\ =
\left. \frac{u^0 z+ u^1 z^2+ u^2 z^3+u^5 z^6}{1-z} 
\exp\left(uz - z + u^2 \frac{z^2}{2} -  \frac{z^2}{2} +
u^3 \frac{z^3}{3} -  \frac{z^3}{3} + u^6 \frac{z^6}{6} -  \frac{z^6}{6}\right) \right|_{u=1} =
\frac{z+z^2+z^3+z^6}{1-z}.$$
The conclusion is that $E[F_6] = 4$ for $n\ge 6$ and there are four fixed points on average.
The general procedure is
$$ G(z, u) = \exp\left(\sum_{d\mid k} 
\left(u^d \frac{z^d}{d} - \frac{z^d}{d}\right) + 
\log \frac{1}{1-z} \right)=
\frac{1}{1-z} \exp\left(\sum_{d\mid k} 
\left(u^d \frac{z^d}{d} - \frac{z^d}{d}\right)\right).$$
Once more continuing as before, we find
$$ \left.\frac{d}{du} G(z, u)\right|_{u=1} =
\left. \frac{\sum_{d\mid k} u^{d-1} z^d}{1-z} 
\exp\left(\sum_{d\mid k} 
\left(u^d \frac{z^d}{d} - \frac{z^d}{d}\right)\right) \right|_{u=1} =
\frac{\sum_{d\mid k} z^d}{1-z}.$$
We have shown that the value of $E[F_k]$ is equal to $\tau(k)$ (the number of divisors of $k$) as soon as $n\ge k.$ It starts out at $1$ for $n=1$ and increases by one every time $n$ hits a divisor of $k$ up to and including $k$ itself.
A: Hint:  what is the probability that $1$ is a fixed point?  Now use the linearity of expectation.
A: In reference to the calculation above I can offer some Maple code that confirms the correctness of the result (answer is $\tau(k)$). This algorithm is very fast and does not iterate over all permutations but only over cycle decompositions, which are treated according to their multiplicity in $S_n$. This makes it possible to compute with values like $n=24,$ which would otherwise be out of reach ($24!$ has $23$ digits).
Here is some sample output:

> seq(v(n, 1), n=1..24);
     1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1

> seq(v(n, 6), n=1..24);
     1, 2, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4

> seq(v(n, 12), n=1..24);
     1, 2, 3, 4, 4, 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6

> seq(v(n, 16), n=1..24);
     1, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 5, 5, 5

> seq(v(n, 24), n=1..24);
     1, 2, 3, 4, 4, 5, 5, 6, 6, 6, 6, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 8

This is the code:

with(group);

pet_cycleind_vrec := 
proc(n)
        local p, s;
        option remember;

        if n=0 then return 1; fi;

        expand(1/n*add(a[l]*pet_cycleind_vrec(n-l), l=1..n));
end;

fp := 
proc(n, p, k)
        local nonfixed, q, j, cyc;

        q := p;

        for j to k-1 do
            q := mulperms(q, p);
        od;

        nonfixed := 0;
        for cyc in q do
            nonfixed := nonfixed + nops(cyc);
        od;

        n - nonfixed;
end;

v :=
proc(n, k)
        option remember;
        local t, el, v, q, cf, d, len, deg, res, cycs;

        res := 0;

        for t in pet_cycleind_vrec(n) do
            el := 1;

            cf := t; cycs := [];
            for v in indets(t) do
                deg := degree(t, v);
                cf := cf/v^deg;

                len := op(1, v);
                if len>1 then
                   for d to deg do
                       cycs := 
                       [op(cycs), [seq(el+q, q=0..(len-1))]];
                       el := el + len;
                   od;
                else
                   el := el + deg;
                fi;

            od;

            res := res + cf*fp(n, cycs, k);
        od;

        res;
end;

