Recall that a derangement is a permutation without fixed points. The problem of computing the number of such derangements on a set of size n is also known as the subfactorial, denoted !n and was first solved by Pierre Raymond de Montmort in 1713. It is now known that

!n = \left\lfloor \frac{n!}{e} \right\rfloor

which is the source of the joke about that the probability that none of n drunk mathematicians get their own coat after a dinner party is 1/e.

A simple way to sample a random derangement is use rejection sampling and sample random permutations until you get a derangement, e.g. using the Fisher-Yates shuffle. So on average, e samplings will be sufficient to get a derangement.

A nicer way to do this, which avoids the rejection sampling, is presented in “Generating Random Derangements” by Martínez, Panholzer and Prodinger. This method has complexity similar to that of the Fisher-Yates shuffle and is simple to implement. Here, I present a slight variation of Martínez et.al’s algorithm with two improvements:

• It avoids a rejection sampling on the index to swap with. I found this a bit annoying since one of the reasons for not using Fisher-Yates is to avoid rejection sampling.
• By using a different data structure for the marked elements, I avoid one of the nested if statements.

The complexity and performance should be about the same as the original algorithm but it is slightly more elegant in my opinion.

The function looks as follows in Rust. It assumes that there is a function, subfactorials, which computes the subfactorials from 0 to n as described above:

pub fn sample_derangement(n: usize) -> Vec<usize> {
    let mut rng = rand::rng();

    let d = subfactorials(n);

    let mut permutation = (0..n).collect::<Vec<usize>>();
    let mut unmarked = (0..n).collect::<Vec<usize>>();

    let mut u = n - 1;
    while u > 0 {
        let i = unmarked.pop().unwrap();
        let j = rng.random_range(0..unmarked.len());
        permutation.as_mut_slice().swap(i, unmarked[j]);
        if rng.random_bool(d[u - 1] as f64 * u as f64 / d[u + 1] as f64) {
            unmarked.remove(j);
            u -= 1;
            if u == 0 {
                break;
            };
        };
        u -= 1;
    }
    permutation
}