note dcturner <p><i>Here is a recent article describing a simpler algorithm for sampling derangements...</i></p> <p>Ah yes, that's nice. I think there is a slight buglet in the given code:</p> <code>\$j = 1 + int rand(\$i-2) while \$mark[\$j];</code> <p>should read</p> <code>\$j = 1 + int rand(\$i-1) while \$mark[\$j];</code> <p>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.</p> <p>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 <i>almost certainly</i> terminates :) On the minus side, it is a bit quadratic (although with low probability). Side-by-side, they run pretty similarly it seems.</p> <code> 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; } </code> 695750 728778