Recurrence for the number of $\sigma \in S_n$ with cycle length at most $r$ I have just learned that the formula is right, but the definition of $c_n^{(r)}$ was wrong. The correct problem is: Prove $$c_{n+1}^{(r)} = \sum_{k=n-r+1}^n n^{\underline{n-k}} c_k^{(r)}$$ where $c_n^{(r)} = |\left\{\sigma \in S_n:\ \text{the length of each cycle of $\sigma$ is at most $r$} \right\}|.$
I will try to prove this with the same approach: classifying by the length of the cycle containing $n+1$.
[edit]
I've just tried this and one details still puzzles me: I assume the cycle that containing $n+1$ to have length $n+1-k$ which leaves $k$ elements for the other cycles. There are $c_k^{(r)}$ ways to make them admissible and I want to accumulate all the cases for $n+1-k$ from 1 to $n$, thus $k$ from $n+1-r$ to $n$. So far, so good. Now we need to consider all possible permutations of the cycle containing $n+1$. But here's my problem: Those should be $(n+1)^{\underline{(n+1)-k}}$ because I have different possibilities for the position of $n+1$ itself. But the formula doesn't seem to account for that.
I have thought about this for quite some time now, but I'm still confused. Can anyone help me answering this?
I also considered $n^{\underline{n-k}}=(n-k)! \binom{n}{n-k}$ which basically gives the same meaning: We choose $n-k$ elements for the cycle containing $n+1$ and consider all possible permutations for them. Clearly, I don't have to choose $n+1$ in the first place, but why is $n+1$ not included in those permutations?
[Old question before correction:]
I want to prove this formula: $$c_{n+1}^{(r)} = \sum_{k=n-r+1}^n n^{\underline{n-k}} c_k^{(r)}$$ where $c_n^{(r)} = |\left\{\sigma \in S_n:\ \sigma^r = id \right\}|$ and I will write $C_n^{(r)} = \left\{\sigma \in S_n:\ \sigma^r = id \right\}$ for the respective set. $n^{\underline{n-k}} = n(n-1)...(n-k+1)$ are the falling factorials and $[n]=\left\{1,...n\right\}.$
I tried simple non-trivial examples and looked at the cycle representation. For $\sigma \in C_n^{(r)}$ the length of every cycle of $\sigma$ has to divide $r$ (of course that condition is not sufficient). But I don't see how that helps so far.
Right now, I'm trying to give some meaning to the right hand side of the recurrence formula and so far I came up with this: Consider the cycle representation of $\sigma \in C_{n+1}^{(r)}$. I can designate the cycle that contains $n+1$. It's length is between $1$ and $r$. Since we have $c_k^{(r)}$ on the right hand side, I want to be left with $k$ elements of $[n+1]$ for the remaining cycles, so the length of the cycle that contains $n+1$ should be $n+1-k$. This also makes sense because it means $1 \leq n+1-k \leq r$ and hence $n+1-r \leq k \leq n$. That's why the indices go from $n+1-r$ to $n$. Now, for $k=n+1-r,...,n$, $c_k^{(r)}$ should be the number of ways to arrange the remaining $k$ elements of $[n+1]$ (those that don't belong to the cycle with $n+1$) into cycles such that $\sigma^r=id$ is possible (necessary condition).
However, at this point I'm stuck because I don't see how the falling factorials fit in, and I don't see how it gives a sufficient condition for $\sigma^r=id$.
 A: I am not sure about your  recurrence, but here is another one that can
be derived.
We are working with the labelled species
$$\mathfrak{P}\left(\sum_{d|r} \mathfrak{C}_d(\mathcal{Z})\right).$$
This yields the generating function
$$G(z) = \exp\left(\sum_{d|r} \frac{z^d}{d}\right).$$
Therefore
$$C_{n+1}^{(r)} =
(n+1)! [z^{n+1}] \exp\left(\sum_{d|r} \frac{z^d}{d}\right).$$
Differentiate to obtain
$$C_{n+1}^{(r)} =
n! [z^n] 
\exp\left(\sum_{d|r} \frac{z^d}{d}\right)
\sum_{d|r} z^{d-1}.$$
What we have here is  the convolution of two generating functions. The
first of $C_n^{(r)}$ and the second of $(d-1)!$. Therefore
$$C_{n+1}^{(r)} =
\sum_{d|r} {n\choose d-1} (d-1)! 
\times(n-d+1)! [z^{n-d+1}] 
\exp\left(\sum_{d|r} \frac{z^d}{d}\right).$$
This finally yields
$$C_{n+1}^{(r)} =
\sum_{d|r} {n\choose d-1} (d-1)! C_{n-d+1}^{(r)}
= \sum_{d|r} \frac{n!}{(n-d+1)!} C_{n-d+1}^{(r)}.$$
Using the notation by Capelli and Toscano as documented
at Wikipedia 
this becomes
$$C_{n+1}^{(r)} =
\sum_{d|r} n^{\underline{d}} C_{n-d+1}^{(r)}.$$
This formula is in fact quite trivial. Simply do the enumeration by classifying according to the cycle that $n+1$ is on.
Addendum.
Observe that when we  multiply two exponential generating functions of
the sequences $\{a_n\}$ and $\{b_n\}$ we get that
$$ A(z) B(z) = \sum_{n\ge 0} a_n \frac{z^n}{n!} 
\sum_{n\ge 0} b_n \frac{z^n}{n!}
= \sum_{n\ge 0} 
\sum_{k=0}^n \frac{1}{k!}\frac{1}{(n-k)!} a_k b_{n-k} z^n\\
= \sum_{n\ge 0} 
\sum_{k=0}^n \frac{n!}{k!(n-k)!} a_k b_{n-k} \frac{z^n}{n!}
= \sum_{n\ge 0} 
\left(\sum_{k=0}^n {n\choose k} a_k b_{n-k}\right)\frac{z^n}{n!}$$
i.e. the  product of  the two generating  functions is  the generating
function of $$\sum_{k=0}^n {n\choose k} a_k b_{n-k}.$$
Addendum Sat Oct 10 CEST 2015. The recurrence above can just as easily be done without the convolution. We have
$$n! [z^n] 
\exp\left(\sum_{d|r} \frac{z^d}{d}\right)
\sum_{d|r} z^{d-1}
= n! \sum_{d|r} [z^{n-d+1}] G(z)
= n! \sum_{d|r} \frac{C_{n-d+1}^{(r)}}{(n-d+1)!}$$
and may continue as before. 
Here the sum is understood to range over those $d$ where $n-(d-1)\ge 0$ as in the following Maple program:

C := proc(n, r)
local m, s, d;
option remember;
    if n = 0 then return 1 end if;
    m := n - 1;
    s := 0;
    for d in numtheory:-divisors(r) do
        if d - 1 <= m then s := s + C(m - d + 1, r)/(m - d + 1)! end if
    end do;
    m!*s
end proc

