Re^3: 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];
}
[download]

blokhead

Comment on Re^3: Random Derangement Of An Array Download Code

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;
}
[download]

[reply]
[d/l]
[select]


Clear questions and runnable code get the best and fastest answer
	PerlMonks