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, 211_F (where the subscript denotes "factorial base") represents decimal 15:

  15 = 2*3! + 1*2! + 1*1!
[download]

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 0_F through 321_F, or 0 through 23 decimal².

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]

All we need now is some function n_to_fb³ 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]

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 100_F from permutation 000_F by swapping position 0 and position 1 (counting from the left). Similarly, one can obtain permutation 200_F from permutation 100_F by swapping position 0 with position 2. The same pattern holds to obtain 300_F from 200_F. The pattern can be expressed succinctly like this:

Rule 1: to obtain permutation b00_F from permutation a00 (where b−a = 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 b00_F we need the permutation for a00_F, 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 *00_F 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]

This is no coincidence. After all, iterating from permutation with rank a00_F to the one with rank a21_F amounts to traversing the lexicographic ordering of the last three letters of a00_F. 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 a00_F. This simple relation between the permutations with ranks a00_F and a21_F enables us to modify our rule above so that it can operate on the immediately preceding permutation:

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

yadda-yadda

(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 xb0_F from permutation xa0_F (where b−a = 1) swap positions 1 and 1+b (counting from the left).

For example, we obtain BCAD (110_F) from BACD (100_F) by swapping positions 1 and 2; and BDAC (120_F) from BCAD (110_F) by swapping positions 1 and 3. (Note that the requirement in Rule 2 stating that b−a 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 xb0_F from permutation xa1_F (where b−a = 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 xyb_F from permutation xya (where b−a = 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]

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]

Final comments

Even though I relied heavily on the notion of lexicographic ordering in my explanation of it, the algorithm itself doesn't do much explicit ordering. It can be said that it iterates in lexicographic order according to the order defined by the input list. If one starts out with the list B, A, D, C, the first few permutations visited will be BADC, BACD, BDAC, BDCA, etc., which is the lexicographic order when B < A < D < C. I see this as a nice bonus feature of the algorithm (in HOP, the problem was stated as simply enumerating all permutations, not in any particular order). But even if one gives it as input the list A, A, A, A, it will dutifully return 24 identical "permutations", AAAA, AAAA, AAAA, etc., and it would be meaningless to call this a lexicographic ordering. Whether this is a feature or a bug is somewhat debatable.

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]

    @arr[ 0..$i-1 ] = reverse @arr[ 0..$i-1 ];
[download]

    if ( $i % 2 ) {
      @arr[ $i-1, 0 ] = @arr[ 0, $i-1 ];
    }
    else {
      @arr[ $i-1, $d-1 ] = @arr[ $d-1, $i-1 ];
    }
[download]

Iterator {...}
[download]

\sum_{k=0}^{m-1}k*k! = m!-1
[download]

(m - 1)*(m - 1)! + (m - 1)! - 1 = m! - 1
[download]

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]

⁴In HOP, Dominus does not mention the Cantor expansion. He refers to the addresses I described above as "patterns", except that his patterns always have an extra 0 tacked on at the end (to "splice" out the last element in a 1-element array).▲