Random Derangement Of An Array

[Limbic~Region]

All,
I need to derange the values of an array randomly (for some definition of random). I could use tye's very nice Derangements iterator and call it a random number of times, but that seems wasteful. I could also use List::Util's shuffle() and continue to swap and elements that remain in their original positions until all are swapped - but that's seems clumsy.

Anyone have a good algorithm for generating a random derangement?

[pc88mxer]

Have a look at this thread: random derangement. Suggestions include:
1. generating a random permutation until you get a derangement (expected number of tries is e), and:
2. using a modification of the basic random shuffle algorithm
It's a Mma mailing list, but I'm sure the algorithms can be translated to perl.

Update: the modification of the basic random shuffle algorithm is to swap position j with a random element from positions j through the end which doesn't violate the derangement condition.

[Corion]

All I'm coming up with are trivial methods:

You can just shift each item by n, where @array > n > 1:

my $shift = rand $#array;
@array = map { $array[ ($_ + $shift) % @array ]} 0..$#array;

Also, reverse is a derangement for even-sized lists, and for odd-sized lists, I guess you can prepend the middle item and reverse the rest:

my $shift = rand $#array;

if (@array % 2 == 0) {
    @array = map { $array[ ($_ + $shift) % @array} reverse 0..$#array;
} else {
    my $middle = int(@array /2);
    my @middle = splice @array, $middle, 1;
    @array = @middle, map { $array[ ($_ + $shift) % @array} reverse 0.
+.$#array;
};

... but all of these derangements are pretty predictable (that is, from the position of one item, you can infer the positions of all other items with two or three guesses)

[Joost]

How 'bout this for an algorithm:

1. shuffle the list using List::Util's shuffle() routine.

2. swap all pairs of elements that are in their original place.

3. figure out what to do when there's an uneven number of elements in their original place (this shouldn't be hard).

[Limbic~Region]

All,
Here is an idea I came up with in #perl (freenode) thanks to a few folks there. Assume an even number of elements:
• 1. Shuffle the indices of the array (copy not actual values)
• 2. Swap adjacent pairs (actually change values)

For instance: 0 - 9 might become 5, 3, 1, 2, 4, 7, 9, 8, 0, 6
Swap index 5 with 3, 1 with 2, 4 with 7, and 9 with 8

[tye]

Shuffle the numbers 1..$#list. Add in the number 0 (just on the end). Then break this list up into 1 or more substrings of at least 2 indices each. Rotate the items at the indices of each substring.

Breaking the list up into pairs (as you did above) is one way to do that if the number of items is even. Perhaps better is to just keep it one list and rotate the whole list, as demonstrated at Re: Derangements iterator.

If you remove the "at least 2 indices each" requirement then rotating gives you permutations, not derangements. Of course, not all derangements (nor permutations) are equally likely when using this technique. But it does select from the subset of derangements that are random rotations uniformly.

Which method is more appropriate depends on the X of your XY problem (which you didn't specify, of course).

[Limbic~Region]

tye,
Which method is more appropriate depends on the X of your XY problem (which you didn't specify, of course).

I am writing a game for a relative to generate cryptograms. Anything not too obvious would be "good enough".

On the other hand, I am interested in more than just "good enough" solutions for the pure pursuit of knowledge.

[BrowserUk]

Try this.

Update: Corrected the rand() calls. It now produces all the possible derangements.

I'm not sure about the statistical performance. Ie. Whether all derangements are equally probably?

#! perl -slw
use strict;
use Data::Dump qw[ pp ]; $Data::Dump::MAX_WIDTH = 1000;

sub derange {
    my $n = shift;
    my @x = 0 .. $n;
    for my $p ( 0 .. $#x ) {
        my $q = int( rand @x );
        $q = int( rand @x ) while $x[ $q ] == $p or $x[ $p ] == $q;
        @x[ $p, $q ] = @x[ $q, $p ];
    }
    return @x;
}

our $N ||= 10;
our $S ||= 1e4;

my @arrange = 0 .. $N;

my %stats;
for( 1.. $S ) {
    my @derange = @arrange[ derange( $#arrange ) ];

    $arrange[ $_ ] == $derange[ $_ ] and warn "Bad at $_\n" for 0 .. $
+#arrange;
    ++$stats{ "@derange" };
#    print "@arrange\n@derange";
}

print "Unique derangements generated: " . keys %stats;
pp \%stats if keys( %stats ) < 50;

my %stats2;
++$stats2{ $_ } for values %stats;

print "Distribution by number of occurances:\n",
    pp \%stats2 if keys( %stats2 ) < 50;

__END__
c:\test>695750 -S=1e3 -N=2
Unique derangements generated: 2
{ "1 2 0" => 651, "2 0 1" => 349 }

Distribution by number of occurances:
{ 349 => 1, 651 => 1 }

c:\test>695750 -S=1e4 -N=3
Unique derangements generated: 9
{ 
    "1 0 3 2" => 1854, "1 2 3 0" => 1744, "1 3 0 2" => 979, 
    "2 0 3 1" => 1028, "2 3 0 1" => 1109, "2 3 1 0" => 1000, 
    "3 0 1 2" => 585, "3 2 0 1" => 1016, "3 2 1 0" => 685 
}

Distribution by number of occurances:
{ 
    585 => 1, 685 => 1, 979 => 1, 1000 => 1, 1016 => 1,
    1028 => 1, 1109 => 1, 1744 => 1, 1854 => 1 
}

c:\test>695750 -S=1e4 -N=4
Unique derangements generated: 44
{ 
    "1 0 3 4 2" => 551, "1 0 4 2 3" => 372, "1 2 0 4 3" => 504, 
    "1 2 3 4 0" => 382, "1 2 4 0 3" => 259, "1 3 0 4 2" => 234, 
    "1 3 4 0 2" => 322, "1 3 4 2 0" => 265, "1 4 0 2 3" => 132, 
    "1 4 3 0 2" => 217, "1 4 3 2 0" => 230, "2 0 1 4 3" => 358, 
    "2 0 3 4 1" => 269, "2 0 4 1 3" => 161, "2 3 0 4 1" => 347, 
    "2 3 1 4 0" => 219, "2 3 4 0 1" => 204, "2 3 4 1 0" => 220, 
    "2 4 0 1 3" => 263, "2 4 1 0 3" => 173, "2 4 3 0 1" => 227, 
    "2 4 3 1 0" => 235, "3 0 1 4 2" => 175, "3 0 4 1 2" => 254, 
    "3 0 4 2 1" => 160, "3 2 0 4 1" => 221, "3 2 1 4 0" => 239, 
    "3 2 4 0 1" => 243, "3 2 4 1 0" => 226, "3 4 0 1 2" => 152, 
    "3 4 0 2 1" => 166, "3 4 1 0 2" => 172, "3 4 1 2 0" => 124, 
    "4 0 1 2 3" => 84, "4 0 3 1 2" => 157, "4 0 3 2 1" => 177, 
    "4 2 0 1 3" => 137, "4 2 1 0 3" => 182, "4 2 3 0 1" => 258, 
    "4 2 3 1 0" => 171, "4 3 0 1 2" => 162, "4 3 0 2 1" => 142, 
    "4 3 1 0 2" => 131, "4 3 1 2 0" => 123 
}

Distribution by number of occurances:
{ 
    84 => 1, 123 => 1, 124 => 1, 131 => 1, 132 => 1, 137 => 1, 
    142 => 1, 152 => 1, 157 => 1, 160 => 1, 161 => 1, 162 => 1, 
    166 => 1, 171 => 1, 172 => 1, 173 => 1, 175 => 1, 177 => 1, 
    182 => 1, 204 => 1, 217 => 1, 219 => 1, 220 => 1, 221 => 1, 
    226 => 1, 227 => 1, 230 => 1, 234 => 1, 235 => 1, 239 => 1, 
    243 => 1, 254 => 1, 258 => 1, 259 => 1, 263 => 1, 265 => 1, 
    269 => 1, 322 => 1, 347 => 1, 358 => 1, 372 => 1, 382 => 1, 
    504 => 1, 551 => 1 
}

c:\test>695750 -S=1e5 -N=5
Unique derangements generated: 265

c:\test>695750 -S=1e5 -N=6
Unique derangements generated: 1854

c:\test>695750 -S=1e5 -N=7
Unique derangements generated: 14544

c:\test>695750 -S=1e6 -N=7
Unique derangements generated: 14833

c:\test>695750 -S=5e6 -N=8
Unique derangements generated: 133496

}

[jethro]

for my $i (1..$#a) {
  my $random= rand( $i);
  swap(\$a[$i], \$a[$random]);
}

Since any position is only swapped once with any of the previous places, there can't be the initial value in any place. But still any value can go to any other position in the array. Per induction I think I can show that the values are equally distributed.

EDIT: Thanks Joost, I was too sloppy, even used { instead of [. If I remember correctly it would work without the \ since @_ in the sub is an alias to the parameters, but this way the side effect is more obvious

[Joost]

I think you meant

 swap(\$a[$i], \$a[$random]);

Which looks like it may work. But randomness always got me confused (and into furious rows with my math teacher, but that's another story).

[blokhead]

I'm skeptical about this one. It's nice, but it can only give factorial(n-1) possible outputs, since that is the number of choices made in the algorithm. But there are many more derangements than that -- the number is essentially factorial(n)/e, which is larger if n>2. So I don't think it gets the entire range of derangements.

For instance, I think one that it wouldn't be able to get is: (4,3,2,1). In the last step, you must swap the 4 into its final position, so the 4 would have to be swapped with the 1 in the first position. But the 1 would definitely no longer be in the first position at that time.

[jethro]

Excellent example. So that algorithm is sadly lacking, and my proof has a hole as well. That happens when you try to proove something in your head instead of really doing the math.

[blokhead]

Here's the approach I wanted to code up, but I've been sufficiently distracted.

The number of derangements is d(n) = (n-1)( d(n-1) + d(n-2) ), and there is a combinatorial proof of this at the wikipedia article. That is, there are (n-1) ways to build an n-derangement out of a (n-1)-derangement, and (n-1) ways to build an n-derangement out of a (n-2)-derangment. Furthermore, these correspond uniquely to all the ways to build an n-derangement.

So here is a way to randomly generate an n-derangement:

• Recursively generate either a (n-1)-derangement or (n-2)-derangement, with probabilities relative to d(n-1) and d(n-2). The base cases are d(1)=0 (no ways to generate a 1-derangement) and d(2)=1 (only 1 choice for a 2-derangement).
• Randomly pick one of the (n-1) ways to generate an n-derangement from the one you have, according to the combinatorial proof.
Unfortunately, I don't have any time to code this up, and the combinatorial proof is not written well. For the case of making an n-derangement out of an (n-1)-derangement, you simply add n to the end, and then swap the last position and a randomly chosen other position. I couldn't quite understand exactly what the wikipedia article was getting at for the other case, though, and it's been too long since I've thought of such things.

Maybe some curious monk can work this into usable code. But this approach is basically recursive: take a random starting derangement and choose a random way to augment it into a larger one. The end result will certainly be randomly distributed, provided you handle the (n-1) and (n-2) cases with the appropriate relative probabilities..

All I managed to get into code was a simple routine to count d(n). It uses a different combinatorial identity that is more amenable to simple computation. It seems like you'd need this to get the relative probabilities for the (n-1) and (n-2) cases to work out:

sub num_d {
    my ($n) = @_;
    return 1-$n if $n < 2;

    my $d = 0;
    $d += (-1)**$_ + ($_-1)*$d for 2 .. $n;

    return $d;
}

[dcturner]

Here's the approach I wanted to code up, but I've been sufficiently distracted.

After much brain-racking I have been unable to come up with a less silly approach than simply encoding the combinatorial proof that d_(n+1) = n (d_n + d_(n-1)) in Perl, which was essentially your idea.

All of the other ideas in this thread so far produce a biased distribution; the 'obvious' modification to the Fisher-Yates shuffle algorithm looked promising but misses some derangements entirely.

It might be straightforward to do it purely iteratively too, but my brain works recursively so my programs do too!

#!/usr/bin/perl

use strict;

sub d_l_rec {
    # Calculates the pair (d_n, d_{n-1})
    #
    # Why calculate the pair? It's O(n) instead of O(2^n)
    # for the 'obvious' recursive algorithm.
    #
    # Why recursive? No need to do it iteratively. 
    my ($n) = @_;
    return 1 if $n < 1;
    return (0, 1) if $n == 1;
    my ($d1, $d2) = d_l_rec($n-1);
    my $d = ($n-1) * ($d1 + $d2);
    return ($d, $d1);
}

sub random_local_derangement {
    # Returns a randomly-chosen local derangement of
    # (0..($n-1)).  A local derangement is a derangement
    # *except* that the last place may be a fixed point.
    #
    # It's 'local' in the sense that, given $n people and a
    # hat of $n tickets you can generate a local derangement
    # by each person in turn pulling tickets out of the hat
    # until they get one that isn't theirs, ensuring that
    # (local to their draw) the permutation looks like it's
    # going to be a derangement. This is how 'Secret Santa'
    # derangements often seem to be organised, and this
    # process sometimes leaves the last person with their
    # own ticket.
    #
    # Unfortunately the 'Secret Santa' approach definitely
    # does not give uniform probabilities to each outcome.
    # This function, on the other hand, does. [At least,
    # it's definitely uniform for $n <= 12 and merely a very
    # good approximation for larger $n.]
    
    my ($n) = @_;
    if ($n == 0) {
        return [];
    }

    my ($i, $threshold);
    # A local $n-derangement is either a full 'total'
    # $n-derangement or else it is a ($n-1)-derangement with
    # a fixed point added at the end. We must choose between
    # these options with appropriate probability weighting.
    # Note that this means that l_k = d_k + d_{k-1} where
    # l_k is the number of local k-derangements and d_k is
    # the number of total k-derangements.
    #
    # Note that d_{12} < 2^31 < d_{13} so we have to be
    # careful of overflow for $n > 13. Fortunately there's a
    # good approximation that can be used.
    if ($n <= 12) {
        # Calculate d_{$n} and d_{$n-1} and a random value
        # $i in [0, d_{$n} + d_{$n-1}) to decide which of
        # the two options to use.
        my ($dn, $dn1) = d_l_rec($n);
        $i = int(rand($dn+$dn1));
        $threshold = $dn1;
    } else {
        # If $n is large then d_$n/d_{$n-1} is very close to
        # $n.  Therefore it'll do to pick a random $i in the
        # range [0, $n] and see if it is 0 or not.
        $i = int(rand($n+1));
        $threshold = 1;

        # I think this is ok, but it just might contain an
        # off-by-one error. The upshot of such an error
        # would be a degree of bias in the results that is
        # going to be hard to detect - you may have to run
        # it literally trillions of times to pick up a
        # statistically significant result.
    }

    if ($i < $threshold) {
        # Case 1 - pick a properly local derangement
        my $d = random_derangement($n-1);
        push @$d, $n-1;
        return $d;
    } else {
        # Case 2 - pick a total derangement
        my $d = random_derangement($n);
        return $d;
    }
}

sub random_derangement {
    # Returns a randomly-chosen (total) derangement of
    # (0..($n-1)), uniformly-chosen amongst all possible
    # derangements.

    my ($n) = @_;
    if ($n == 0) {
        return [];
    }

    # There are (n-1) l_{n-1} of them, so pick a (uniformly)
    # random local ($n-1)-derangement and a random $m in the
    # range [0, $n-1).
    my $ld = random_local_derangement($n-1);
    my $m = int(rand($n-1));

    # If L_k is the set of all local k-derangements and D_k
    # is the set of all total k-derangements then the code
    # below encodes the proof that (n-1) l_{n-1} = d_n in a
    # bijection between [0, $n-1) x L_{n-1} and D_{n}.
    #
    # Since the pair ($m, $ld) are chosen uniformly, this
    # shows that the resulting derangement is also uniformly
    # chosen.

    if ($n-2 == $ld->[$n-2]) {
        # $ld is properly local. Therefore the desired
        # derangement swaps the $m'th and last places and
        # uses $ld to derange the other places.

        my $j = $n-1;
        while ($j--) {
            my $k = $j < $m ? $j : $j-1;
            $ld->[$j] 
                = $ld->[$k] < $m ? $ld->[$k] : $ld->[$k]+1;
        }
        $ld->[$n-1] = $m;
        $ld->[$m] = $n-1;
        return $ld;
    } else {
        # $ld is total. Therefore put the $m'th entry at the
        # end and put $n-1 in the $m'th place.

        $ld->[$n-1] = $ld->[$m];
        $ld->[$m] = $n-1;
        return $ld;
    }
}

sub check_derangement {
    # Check that we have really generated a derangement
    my ($n, $d) = @_;

    my $s = join ', ', @$d;

    die "Wrong length: $s ($n)" unless ($n == @$d);
    for (my $i = 0; $i < $n; $i++) {
        die "Not a derangement: $s" if ($i == $d->[$i]);
        die "Illegal value: $d->[$i] in $s" 
            if ($d->[$i] < 0 || $d->[$i] >= $n);
    }

    eval {
        my @check_unique = sort { $a <=> $b || undef } @$d;
    };
    die "Uniqueness check failed: $s" if $@;
}

my $n = $ARGV[0];

my %f = ();
my $c = 0;

for (my $i = 0; $i < 1e6; $i++) {
    my $d = random_derangement($n);
    check_derangement($n, $d);
    my $s = join ',', @$d;
    $f{$s} += 1;
    $c += 1;
}

for my $key (sort {$f{$a} <=> $f{$b}} keys %f) {
    printf "%s: %0.2f%% (expected %+0.2f%%)\n", $key, 100.0*($f{$key}/
+$c),
        100.0*(($f{$key}/$c)-(1.0/(scalar keys %f)));
}

printf "Total %d (%d runs)\n", (scalar keys %f), $c;

This produces output like the following

$ ./derangements.pl 4
2,3,0,1: 11.08% (expected -0.03%)
1,0,3,2: 11.09% (expected -0.02%)
2,3,1,0: 11.10% (expected -0.01%)
1,3,0,2: 11.11% (expected -0.01%)
2,0,3,1: 11.11% (expected -0.00%)
1,2,3,0: 11.12% (expected +0.01%)
3,0,1,2: 11.12% (expected +0.01%)
3,2,1,0: 11.13% (expected +0.02%)
3,2,0,1: 11.14% (expected +0.03%)
Total 9 (1000000 runs)

Clearly that's not far off uniform!

[blokhead]

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

[jdporter]

sub random_derangement
{
    my @i = shuffle( 0 .. $#_ );
    my @r;
    @r[ @i ] = @_[ @i[1..$#i], $i[0] ];
    @r
}

If you want it to derange the list in place, it's even simpler:

sub derange
{
    my @i = shuffle( 0 .. $#_ );
    @_[ @i ] = @_[ @i[1..$#i], $i[0] ];
}

Update: I just realized that this is essentially the idea tye suggested in Re^2: Random Derangement Of An Array (rotate). The difference is that I rotate the entire list of (shuffled) indices. I don't believe breaking up the list into smaller chunks is necessary, and, depending on how that was done, could actually hurt the randomization a bit.

Update Again: Upon further reflection, I believe tye is right. In my algorithm above, some derangements are impossible. I had originally had the same idea — subdividing the set of indices — before I (incorrectly, I now believe) made the above oversimplification. The problem is that I don't know how to partition the set fairly. I did come up with an algorithm, but my intuition says it's not exactly fair either. Here it is; perhaps someone can say how fair it is:

sub random_derangement
{
    my @i = shuffle( 0 .. $#_ );
    my @j = @i;

    my @part; # aoa; will contain the partitions
    # distribute the elements of @j across the existing partitions ran
+domly.
    # but if there are N partitions, there is a 1/(N+1) chance that th
+e element
    # will be distributed to a new (N+1)th partition.
    # if the number of existing partitions containing only a single el
+ement
    # is equal to the number of remaining elements of @j, we can't cho
+ose 
    # just any partition; we have to distribute the remaining elements
+ of
    # @j to each of those existing partitions having a single element.
    # furthermore, we take a precaution against getting into a situati
+on
    # where we'll have more "singleton" partitions than we have remain
+ing
    # elements in @j
    while (@j)
    {
        my @d = grep { @$_ == 1 } @part;
        @d = shuffle( @d );
        if ( @d == @j )
        {
            push @{  $d[0]  }, shift @j;
        }
        elsif ( @d+1 >= @j )
        {
            push @{ $part[ rand( @part ) ] }, shift @j;
        }
        else
        {
            push @{ ($part[ rand( 1 + @part ) ]||=[]) }, shift @j;
        }
    }
    # now do the rotations, re-using @i and @j.
    @i=();
    foreach my $part ( @part )
    {
        push @i, @$part;
        push @j, @{$part}[ 1 .. $#{$part}, 0 ];
    }

    my @r;
    @r[ @i ] = @_[ @j ];
    @r
}

[jacques]

sub fisher_yates_shuffle {
        my $deck = shift;  # $deck is a reference to an array
        my $i = @$deck;
        while ($i--) {
            my $j = int rand ($i+1);
            @$deck[$i,$j] = @$deck[$j,$i];
        }
    }

my @list = (1..5);
fisher_yates_shuffle( \@list );
print join ', ',@list;

[Limbic~Region]

jacques,
Is this a modified Fisher-Yates shuffle? I am too tired to tell. If it isn't, then it won't work because a shuffle is not guaranteed to generate a derangement. If it is - please explain the modification.

[Anonymous Monk]

http://search.cpan.org/~allenday/Math-Combinatorics-0.09/lib/Math/Combinatorics.pm#derange()?

[BrowserUk]

Shame that, according to its own documentation, it doesn't work!

 Usage   : my @deranges = derange(@n);
 Function: implements !n, a derangement of n items in which none of th
+e
           items appear in their originally ordered place.
 Example : my @n = qw(a b c);
           my @d = derange(@n);
           print join "\n", map { join " ", @$_ } @d;
           # prints:
           # a c b  ### Not a derangment; a in it's original position
           # b a c  ### Not a derangment; c in it's original position
           # b c a  
           # c a b 
           # c b a  ### Not a derangment; b in it's original position

[ysth]

I could also use List::Util's shuffle() and continue to swap and elements that remain in their original positions until all are swapped - but that's seems clumsy.

I have the feeling that would introduce some bias - which may not be important for your problem.

[Anonymous Monk]

The generation of random numbers is far too important to be left to chance.