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$.
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