Perl Monk, Perl Meditation PerlMonks

### Re^3: Random Derangement Of An Array

 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];

\$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;
}

Create A New User
Node Status?
node history
Node Type: note [id://728778]
help
Chatterbox?
 [zentara]: I'm listening to 20,000 Tibetan monks chant Om Mani Padme Hum [zentara]: I've decided the first words I will utter when I'm born into that other world..... Hari Om :-)

How do I use this? | Other CB clients
Other Users?
Others wandering the Monastery: (7)
As of 2017-05-24 21:47 GMT
Sections?
Information?
Find Nodes?
Leftovers?
Voting Booth?
My favorite model of computation is ...

Results (186 votes). Check out past polls.