Beefy Boxes and Bandwidth Generously Provided by pair Networks
Perl: the Markov chain saw

Re^3: Random Derangement Of An Array

by blokhead (Monsignor)
on Dec 07, 2008 at 20:38 UTC ( #728778=note: print w/replies, xml ) Need Help??

in reply to Re^2: Random Derangement Of An Array
in thread Random Derangement Of An Array

Here is a recent article describing a simpler algorithm for sampling derangements. I also found slides for the presentation. Since the paper is so recent, I guess this means that a small modification of Fisher-Yates is unlikely to generate derangements, since someone would have already come up with it by now. Still, their algorithm is in-place and has better expected running time than retrying Fisher-Yates until you get a derangement.

Here is a Perl implementation I whipped up. It is slightly odd because I followed their lead and used array indexing from 1.

sub rand_derangement { my $n = shift; return if $n == 1; ## no derangements of size 1 ## precompute $D[n] == number of derangements of size n my @D = (1,0); push @D, $#D * ($D[-1] + $D[-2]) while $#D < $n; my @A = (undef, 1 .. $n); my @mark = (1, (0) x $n); my ($i, $u) = ($n, $n); while ($u > 1) { if (! $mark[$i]) { my $j = 0; $j = 1 + int rand($i-2) while $mark[$j]; @A[$i,$j] = @A[$j,$i]; if ( rand(1) < ($u-1) * $D[$u-2] / $D[$u] ) { $mark[$j] = 1; $u--; } $u--; } $i--; } return @A[1..$n]; }


Replies are listed 'Best First'.
Re^4: Random Derangement Of An Array
by dcturner (Initiate) on Dec 08, 2008 at 19:11 UTC

    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($i-2) while $mark[$j];

    should read

    $j = 1 + int rand($i-1) 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 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; }

Log In?

What's my password?
Create A New User
Node Status?
node history
Node Type: note [id://728778]
[AnomalousMonk]: choroba: It even has Rama's (relatively) high rate of spin, although I don't know if it's axial spin.
[LanX]: Rama spins?
LanX probably in my stomach
erix considers showing this to choroba , thinks better of it

How do I use this? | Other CB clients
Other Users?
Others wandering the Monastery: (10)
As of 2017-11-22 16:59 GMT
Find Nodes?
    Voting Booth?
    In order to be able to say "I know Perl", you must have:

    Results (327 votes). Check out past polls.