Derivation of the Partial Derangement (Rencontres numbers) formula

I'm looking for the method by which the partial derangement formula $D_{n,k}$ was derived. I can determine the values for small values of N empirically, but how the general case formula arose still alludes me.

Any links/books or an explanation will be appreciated.

The formula is:

$D_{n,k} = {n \choose k}!(n-k)$

Links:

1. Mathworld

Here is a very general solution. There is a fundamental formula in combinatorics called the exponential formula, and one statement of it is as follows. Given a finite group $G$ acting on a set $X$, its cycle index polynomial is given by

$$Z_G = \frac{1}{|G|} \sum_{g \in G} z_1^{c_1(g)} z_2^{c_2(g)} ... $$

where $c_i(g)$ is the number of cycles of length $i$ in the action of $g$ on $X$. In particular, the notation $Z_{S_n}$ will denote the cycle index polynomial of $S_n$ acting on an $n$-element set in the usual way; it is a generating function encoding the relative proportions of different cycle types of permutations.

The exponential formula then states that

$$\sum_{n \ge 0} Z_{S_n} t^n = \exp \left( z_1 t + \frac{z_2 t^2}{2} + \frac{z_3 t^3}{3} + ... \right).$$

In my opinion this is one of the most beautiful formulas in mathematics and a major reason I became interested in combinatorics was because I stumbled upon this formula while solving a Putnam problem (which is described in the blog post I linked to above).

How does it apply to this problem? Set $z_2 = z_3 = ... = 1$ and $z_1 = z$. Then the LHS of the exponential formula is a generating function which counts permutations according to how many fixed points ($1$-cycles) they have. In other words,

$$Z_{S_n}(z, 1, 1, ...) = \frac{1}{n!} \sum_{g \in S_n} z^{c_1(g)} = \frac{1}{n!} \sum_{k=0}^n D_{n,k} z^k.$$

The RHS of the exponential formula, on the other hand, is

$$\exp \left( zt + \log \frac{1}{1-t} - t \right) = \frac{e^{-t}}{1 - t} e^{zt}.$$

So we obtain the beautifully concise formula

$$\sum_{n \ge 0} \frac{t^n}{n!} \sum_{k=0}^n D_{n,k} z^k = \frac{e^{-t}}{1 - t} e^{zt}.$$

The coefficients of $\frac{e^{-t}}{1 - t}$ are obtained by setting $z = 0$; they give the usual derangement numbers, e.g. the number of permutations of $n$ elements with no fixed points, and this can also be seen directly from the generating function since

$$\frac{e^{-t}}{1 - t} = \sum_{n \ge 0} \left( \sum_{k=0}^n \frac{(-1)^k}{k!} \right) t^n$$

which is equivalent to the formula $D_{n,0} = n! \sum_{k=0}^n \frac{(-1)^k}{k!} \sim \frac{n!}{e}$. (In fact you can read this asymptotic directly from the generating function.) The above then gives

$$D_{n,k} = {n \choose n-k} D_{n-k,0} = \frac{n!}{k!} \sum_{i=0}^{n-k} \frac{(-1)^i}{i!}.$$

Of course, there is a much more direct proof of this: observe that specifying a permutation of $n$ elements with $k$ fixed points is equivalent to specifying the $n-k$ elements which are not fixed points, then specifying a fixed-point-free permutation of these. This immediately gives $D_{n,k} = {n \choose n-k} D_{n-k,0}$, so it suffices to compute $D_{n,0}$, and this can be done by the standard inclusion-exclusion argument.

(In the interest of completeness, the standard inclusion-exclusion argument is as follows: first we start with all $n!$ permutations. Then we subtract the ones which fix $1$, and the ones which fix $2$, etc., so we subtract $n \cdot (n-1)!$. But this is overcounting: we need to add back the ones which fix both $1$ and $2$, or more generally both $i$ and $j$ for distinct $i, j$, so we add back ${n \choose 2} \cdot (n-2)!$. But this is overcounting: we need to subtract the ones which fix any triplet... and so forth. This gives each term of the formula $n! \sum_{k=0}^n \frac{(-1)^k}{k!}$ one-by-one.)

My point in presenting the generating function argument is not that it is any easier in this case but that it generalizes to far more complicated problems in a way which minimizes mental effort: for example you can use it to count permutations by how many $2$-cycles they have, or by $c_3(g) + 17 c_5(g)$, or whatever, and the generating function is also a convenient way to organize the computation of the expected value and variance of various permutation statistics; see, for example, this math.SE answer.

The answer by Qiaochu Yuan is excellent, but perhaps admits some simplification. Note that the combinatorial species $\mathcal{Q}$ of permutations by cycles with the fixed points marked is given by $$\mathcal{Q} = \mathfrak{P}(\mathcal{U}\mathfrak{C_1}(\mathcal{Z}) + \mathfrak{C_2}(\mathcal{Z}) + \mathfrak{C_3}(\mathcal{Z}) +\cdots)$$ which immediately yields the generating function $$Q(z, u) = \exp\left(uz + \frac{z^2}{2} + \frac{z^3}{3} + \cdots\right) = \exp\left(uz -z + \log\frac{1}{1-z}\right) = \frac{1}{1-z} \exp\left(uz-z\right).$$ The coefficients from this generating function produce the rencontres numbers as in $$D_{n,k} = n! [z^n] [u^k] Q(z, u).$$ Actually doing the extraction we obtain $$D_{n,k} = n! [z^n] \frac{1}{1-z} \exp(-z) \frac{z^k}{k!} = \frac{n!}{k!} [z^{n-k}] \frac{1}{1-z} \exp(-z) = \frac{n!}{k!} \sum_{m=0}^{n-k} \frac{(-1)^m}{m!}.$$ Maybe this can contribute to an understanding of the essentials of this computation.