On pp. 134-135 of HOP there's a nifty algorithm I had never seen before for iterating over all permutations of a list (though I've seen some that look similar). It looked to me like the full explanation for this algorithm ended up on the editing room floor (the brief explanation given in the book is very unclear), plus I was not able to readily find one online, so I figure that a little clarification of this algorithm may be a good thing to have. Also, in the process of understanding this algorithm one comes across some cool math.

To pique your interest, here's the algorithm in its original form1:

sub permute {
my @items = @_;
my $n = 0; return sub {$n++, return @items if $n==0; my$i;
my $p =$n;
for ($i=1;$i<=@items && $p%$i==0; $i++) {$p /= $i; } my$d = $p %$i;
my $j = @items -$i;
return if $j < 0; @items[$j+1..$#items] = reverse @items[$j+1..$#items]; @items[$j,$j+$d] = @items[$j+$d,$j];$n++;
return @items;
}
}
[download]
The function permute just returns an iterator, which is where the action is. The first time around the iterator simply returns the input list, but after that things get more interesting.

(It's a nice puzzle to figure out why this algorithm works. If you want to figure it out by yourself, stop reading here.)

#### Lexicographic order, revisited

We're all familiar with lexicographic sorting; in fact that's what perl's sort does by default. But bear with me as I describe in a somewhat roundabout way the procedure for sorting lexicographically all the permutations ABCD, ABDC, ACBD, etc. of the letters A, B, C, and D.

First we group all the permutations into 4 groups, according to their first letter, and we number these 4 groups, starting with 0, according to the alphabetic order of this first letter. Hence, group 0 is the one for the permutations beginning with A, etc. Now, for each of these groups we perform the same group-and-number procedure, using the second letter as the basis for the grouping. In this case the number of groups is 3, since, within each of the top-level groups, only three letters are available for the second position. Then we repeat this process one more and last time: grouping according to the third letter (only two are possible at this level) and numbering the resulting groups according to the alphabetic order in the third position.

Now each permutation belongs to three nested (and numbered) groups. We can use the numbers of these groups to assign a unique "address" to each permutation. For example, the permutation CADB has address 201, because the C's are group 2 at the top level, the CA's are group 0 within the C's, and the CAD's (which contain only one permutation, CADB) is group 1 within the CA's. By this addressing scheme, the 24 permutations get addresses ranging from 000 for ABCD to 321 for DCBA.

#### The Cantor expansion

There is a straightforward one-to-one mapping between this addressing scheme and the natural numbers; it's called the Cantor expansion, or factorial base expansion. For example, 211F (where the subscript denotes "factorial base") represents decimal 15:

  15 = 2*3! + 1*2! + 1*1!
[download]
The principle of the representation is similar to that of the familiar representation of integers as powers of a fixed base, excepts that the digits in the factorial-base reprsentation correspond to coefficients of factorials.

If we interpret the addresses we derived earlier (from the lexicographic sorting) as base-F numbers, then the lexicographic order for the permutations corresponds to numeric order of their addresses, which range from 0F through 321F, or 0 through 23 decimal2.

But more importantly, the addresses encode a simple procedure for mapping a natural number to a given permutation. Here it is:

sub fb_to_perm {
my @arr = @{ shift() };
my $fb = shift; return unless defined$fb;
my @ret;
push @ret, splice @arr, $_, 1 for @$fb;
push @ret, @arr;
return @ret;
}
[download]
The first argument to fb_to_perm is a reference to the array we want to permute, and the second one points to an array whose elements are the factorial-base representation digits of an integer between 0 and m!−1, where m is the length of the list to be permuted. To construct the permutation correspondiing to 201F, fb_to_perm sequentially splices out 2, 0, and 1 from the array (i.e. B, D, C ) and appends what's left, the single letter B (which is completely determined as the fourth position at this point).

All we need now is some function n_to_fb3 for computing the factorial-base representation of a natural number. Then, to iterate over all the permutations of a list, we simply do something like this:

my $n = 0; { my$fb = n_to_fb( $n, scalar @my_list ); last unless$fb;
my @perm = fb_to_perm( \@my_list, $fb ); print "@perm\n"; ++$n;
redo;
}
[download]
In HOP Dominus presents a very similar strategy on pp. 130-134, including an algorithm (pattern_to_permutation4 on p. 131) which is very similar to the one used by fb_to_perm. Ultimately, however, he rejects it because it requires too much splicing, and proposes the algorithm that motivated this meditation as a more efficient (or at least less splice-heavy) alternative.

How does this alternative algorithm manage to avoid all the splicing? The trick is that, instead of building the permutation from scratch every time, it computes it as an in-place transformation of the previous one, using swaps and flips (which preserve the length of the array) instead of splice. The trickiest part is figuring out where to flip and what to swap. It turns out that the Cantor-expansion representation plays a central role in the splice-less variant too.

#### Cutting out splice

Below I list all the permutations of the letters A, B, C, D along with their rank (base-F), in lexicographic order.

A B C D 000_F   B A C D 100_F   C A B D 200_F   D A B C 300_F
A B D C 001_F   B A D C 101_F   C A D B 201_F   D A C B 301_F
A C B D 010_F   B C A D 110_F   C B A D 210_F   D B A C 310_F
A C D B 011_F   B C D A 111_F   C B D A 211_F   D B C A 311_F
A D B C 020_F   B D A C 120_F   C D A B 220_F   D C A B 320_F
A D C B 021_F   B D C A 121_F   C D B A 221_F   D C B A 321_F
[download]

Let's focus for a minute on the permutations whose ranks end in 00, and which are listed along the top row above. Notice that one can obtain permutation 100F from permutation 000F by swapping position 0 and position 1 (counting from the left). Similarly, one can obtain permutation 200F from permutation 100F by swapping position 0 with position 2. The same pattern holds to obtain 300F from 200F. The pattern can be expressed succinctly like this:

Rule 1: to obtain permutation b00F from permutation a00 (where ba = 1) swap positions 0 and b (counting from the left).

This tells us that if we can come up with an algorithm for iterating over the permutations of 3 letters in lexicographic order, we can use Rule 1 above to extend this algorithm to iterate over the permutations of 4 letters. In other words, we first use the 3-letter algorithm for iterating over all the permutations of BCD, and append A at the beginning of each. Then we use Rule 1 to move from move from the first A-permutation (ABCD) to the first B-permutation (BACD), and invoke the 3-letter algorithm again. Etc.

There's a problem with Rule 1, though. To get the permutation for b00F we need the permutation for a00F, which was visited several iterations before (6 to be exact). This is very inconvenient, because it requires keeping track of more permutations than just the previous one. But take a look at the permutations immediately preceding the *00F permutations:

A B C D 000_F

A D C B 021_F
B A C D 100_F

B D C A 121_F
C A B D 200_F

C D B A 221_F
D A B C 300_F
[download]
Notice that we can always recover the permutation with index a00F from the permutation immediately preceding the one with index b00F, which has the representation a21F, by flipping (i.e. reversing) the last three letters in it. For example, by flipping the last three letters of the permutation with rank 121F (BDCA) we get BACD, which is the permutation that is required by Rule 1 to generate the permutation with rank 200F, namely the one with rank 100F.

This is no coincidence. After all, iterating from permutation with rank a00F to the one with rank a21F amounts to traversing the lexicographic ordering of the last three letters of a00F. Therefore, the last permutation in this sub-iteration has to be the one in which the last three letters are in reverse order relative to their order in a00F. This simple relation between the permutations with ranks a00F and a21F enables us to modify our rule above so that it can operate on the immediately preceding permutation:

Rule 1': to obtain permutation b00F from permutation a21 (where ba = 1), reverse the order of the last three letters, and swap positions 0 and b (counting from the left).

In contrast to Rule 1, Rule 1' is completely "local". It tells us how we get certain permutations from the immediately preceding ones, and serves as the model for all other transformation rules we will need.

(In this section we work out the counterparts of Rule 1' for the two remaining cases. If you got it already, skip to the next section, The Mother of All Rules.)

The same idea works at the next level too. Consider the permutations whose indices end in 0.

A B C D 000_F
A C B D 010_F
A D B C 020_F

B A C D 100_F
B C A D 110_F
B D A C 120_F

C A B D 200_F
C B A D 210_F
C D A B 220_F

D A B C 300_F
D B A C 310_F
D C A B 320_F
[download]

Rule 2: to obtain permutation xb0F from permutation xa0F (where ba = 1) swap positions 1 and 1+b (counting from the left).

For example, we obtain BCAD (110F) from BACD (100F) by swapping positions 1 and 2; and BDAC (120F) from BCAD (110F) by swapping positions 1 and 3. (Note that the requirement in Rule 2 stating that ba must equal 1 means that b cannot be 0.)

Therefore, reasoning as before, if we can come up with an algorithm for iterating over the permutations of 2 letters in lexicographic order, we can use Rule 2 above to extend this algorithm to iterate over the permutations of 3 letters (and therefore, by the argument made earlier, we then also have an algorithm for iterating over all permutations of 4 letters).

As with Rule 1, Rule 2 inconveniently requires a permutation other than the previous one, but the fix is just like the fix of Rule 1, except that only the last to letters need to be reversed:

Rule 2': to obtain permutation xb0F from permutation xa1F (where ba = 1), reverse the order of the last two letters, and swap positions 1 and 1+b (counting from the left).

Not surprisingly a similar rule works at the lowest level too.

Rule 3: to obtain permutation xybF from permutation xya (where ba = 1), reverse the last one letter, and swap positions 2 and 2+b (counting from the left).

I included the no-op "reverse the last one letter" because I intend to turn the three rules into one single rule.

#### The Mother Of All Rules

OK, we can now generalize the rule for flipping and swapping.

Flip and Swap Rule: to obtain permutation from the previous one, first determine the number i−1 of trailing zeros and the last non-zero digit d in the factorial-base representation of its lexicographic order rank. Then reverse the order of the last i letters of the previous permutation. Equivalently, reverse the order of the letters at positions j+1 through m−1, where m is the length of the list being permuted and j is defined as m−i. Finally, swap the letters at positions j and j+d.

The most interesting task implied by this rule is determining the numbers i and d corresponding to a given natural number n. The following subroutine takes care of this:

sub offsets {
my $p =$_[ 0 ];
my $i = 2; my$d;
$p /=$i++ until $d =$p % $i; return ($i, $d ); } [download] This little sub is the heart of the whole permutation enumeration algorithm. It takes an integer as input. It divides it by increasingly large divisors (in$i), starting with 2, and incrementing the divisor by 1 each time, as long as this division has 0 remainder. This amounts to counting the number of zeros in the factorial base expansion of $_[ 0 ]. (Actually, it counts 1 more than the number of zeros, which is convenient, because that's what our unified rule requires.) And the first non-zero remainder, which is captured in$d, is the last non-zero digit in this factorial base expansion.

#### The home stretch

Now we can re-construct the original algorithm. The subroutine offsets replaces the for loop inside the original algorithm. It incorporates two very minor optimizations over the version given in HOP. One is that $i starts out with value 2, since starting with value 1 simply wastes an iteration without any effect ($p % 1 is always 0). The other small improvement is that the loop in offsets doesn't bother to check at every iteration whether $i is greater than the length of the original array. Since this events signals the end of the iteration, it is handled by the calling program. In fact, offsets doesn't even know about the original array; its only argument is an integer. sub perms { my @arr = @_; my$m = @arr;
my $n = 0; { print "@arr\n"; ++$n;
my ( $i,$d ) = offsets( $n ); last if$i > @arr;
my $j = @arr -$i;
@arr[ $j+1..$#arr ] = reverse @arr[ $j+1..$#arr ];
@arr[ $j,$j+$d ] = @arr[$j+$d,$j ];
redo;
}
}
[download]

This algorithm is a special case of a large class of algorithms that Don Knuth calls the general permutation generator. Something equivalent to the offsets subroutine above is at the heart of all these generators. The various generators differ in how they order the permutations, which is implied in how the offsets i and d are mapped to a transformation of the previous permutation. For example, one can obtain a different generator, if, in perms, one replaces the definition of $j and the "flip-and-swap" lines  my$j = @arr - $i; @arr[$j+1..$#arr ] = reverse @arr[$j+1..$#arr ]; @arr[$j, $j+$d ] = @arr[ $j+$d, $j ]; [download] with the single flip  @arr[ 0..$i-1 ] = reverse @arr[ 0..$i-1 ]; [download] ...though, as one would expect, this generator does not produce the permutations in lexicographic order. This scheme doesn't even need$d. The scheme Knuth thinks may be most efficient is one proposed by B. R. Heap in 1963, which can be obtained by replacing the same three lines of perms as before with
    if ( $i % 2 ) { @arr[$i-1, 0 ] = @arr[ 0, $i-1 ]; } else { @arr[$i-1, $d-1 ] = @arr[$d-1, $i-1 ]; } [download] Now, instead of a flip and a swap we have only a test for oddness and a swap. The basis for this general framework is a representation scheme for permutation groups called Sims tables. 1Well, almost; in the original, what permute returns is something that looks like this Iterator {...} [download] but as it turns out, Iterator, which is defined on p. 123, is just an identity function for code refs, so I replaced it with the more transparent sub keyword. 2In the general case, the addresses range from 0...00F to q...4321F, where q is 1 less than the length of the list being permuted. This suggests a nice factorial identity: \sum_{k=0}^{m-1}k*k! = m!-1 [download] for all integers m > 0. The inductive proof of this identity is a one-liner: (m - 1)*(m - 1)! + (m - 1)! - 1 = m! - 1 [download] 3Here's a version of n_to_fb; it is almost identical to HOP's n_to_pat (p. 133): sub n_to_fb { my ($n, $m ) = @_; my @fb; for my$i ( 2 .. $m ) { my$r = $n %$i;
unshift @fb, $r;$n = int( $n/$i );
++$i; } return$n ? undef : \@fb;
[download]