Here is a recent article describing a simpler algorithm for sampling derangements...

Ah yes, that's nice. I think there is a slight buglet in the given code:

should read

should it not? Also it could hit an overflow bug for $n > 12, because d_n gets quite large quite quickly. Even with 64-bit ints you overflow for $n > 20. The same approximation works for your method as for mine: for $n > 12 the value (n-1)d_{n-2} / d_n is very close to 1/n.

The iterative version of the recursive program that I wrote above is as follows. It's a bit convoluted because the recursion d_n = (n-1) (d_{n-1} + d_{n-2}) is second-order, so you have to work a bit harder than for a first-order recursion. It's also in-place and requires only a constant amount of auxiliary storage. Plus it certainly terminates, whereas the rejection method merely almost certainly terminates :) On the minus side, it is a bit quadratic (although with low probability). Side-by-side, they run pretty similarly it seems.

sub random_derangement {
    my ($n) = @_;

    my @d = (1, 0, 1, 2, 9, 44, 265, 1854, 14833, 133496,
        1334961, 14684570, 176214841);

    my @t;
    my $i = $n;
    while ($i--) {
        my $m = int(rand($i));
        $t[$i] = $m;
        if ($i <= 12) {
            $i -= 1 if (int(rand($d[$i]+$d[$i-1])) < $d[$i-1]);
        } else {
            $i -= 1 if (int(rand($i+1)) == 0);
        }
    }
    
    for ($i = 0; $i < $n; $i++) {
        if (defined($t[$i])) {
            my $m = $t[$i];
            $t[$i] = $t[$m];
            $t[$m] = $i;
        } else {
            my $j = $i+1;
            my $m = $t[$i+1];
            while ($j--) {
                my $k = $j < $m ? $j : $j-1;
                $t[$j] 
                    = $t[$k] < $m ? $t[$k] : $t[$k]+1;
            }
            $t[$i+1] = $m;
            $t[$m] = $i+1;
            $i += 1;
        }
    }

    return \@t;
}
[download]