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:
$j = 1 + int rand($i2) while $mark[$j];
should read
$j = 1 + int rand($i1) while $mark[$j];
should it not? Also it could hit an overflow bug for $n > 12, because d_n gets quite large quite quickly. Even with 64bit ints you overflow for $n > 20. The same approximation works for your method as for mine: for $n > 12 the value (n1)d_{n2} / 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 = (n1) (d_{n1} + d_{n2}) is secondorder, so you have to work a bit harder than for a firstorder recursion. It's also inplace 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). Sidebyside, 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[$i1])) < $d[$i1]);
} 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 : $j1;
$t[$j]
= $t[$k] < $m ? $t[$k] : $t[$k]+1;
}
$t[$i+1] = $m;
$t[$m] = $i+1;
$i += 1;
}
}
return \@t;
}
