In an earlier post I described a way to sample a uniformly random derangement, a permutation with no fixed points, in linear time and without any big-integer arithmetic. Since then I have created samplers for a couple of similar objects, like involutions, which are permutations equal to their own inverse, and perfect matchings, which are involutions with no fixed points. Writing the third one I noticed that they are all the same algorithm.

The thing they have in common is the cycle structure. A derangement is a permutation whose cycles all have length \geq 2, an involution has cycles of length \leq 2 and a matching has cycles of length exactly 2. In every case we pick a set S of permitted cycle lengths and ask for a uniformly random permutation built from those alone.

Counting such permutations is a single recurrence. Let a_n be the number of permutations of n elements whose cycle lengths all lie in S. The largest element must sit in a cycle of some length k \in S. There are (n-1)(n-2)\cdots(n-k+1) ways to line up its k-1 cycle-mates in order, and the remaining elements form a smaller instance of the same problem, so

a_n = \sum_{k \in S} (n-1)(n-2)\cdots(n-k+1)\, a_{n-k}.

To sample a uniformly random permutation of the given type, repeatedly take the largest remaining element and give it a cycle of length k with probability proportional to the k-th term above, then draw its k-1 partners uniformly at random. Each valid permutation then turns up with probability exactly \frac{1}{a_n}.

One problem is that a_n grows like n!, so forming those weights directly overflows almost immediately. We sidestep this by never building the numbers at all. The probabilities are computed iteratively, so no big integers are needed however large n gets.

The code can be found in the rand_combinatorics crate.