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))$?
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.
Hint: what is the probability that $1$ is a fixed point? Now use the linearity of expectation.
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;