python - Heap's algorithm permutation generator

Question

I need to iterate over permutations of a tuple of integers. The order has to be generated by swapping a pair of elements at each step.

I found the Wikipedia article (http://en.wikipedia.org/wiki/Heap%27s_algorithm) for Heap's algorithm, which is supposed to do this. The pseudocode is:

procedure generate(n : integer, A : array of any):
    if n = 1 then
        output(A)
    else
        for i := 1; i ≤ n; i += 1 do
            generate(n - 1, A)
            if n is odd then
                j ← 1
            else
                j ← i
            swap(A[j], A[n])

I tried to write a generator for this in python:

def heap_perm(A):
    n = len(A)
    Alist = [el for el in A]
    for hp in _heap_perm_(n, Alist):
        yield hp


def _heap_perm_(n, A):
    if n == 1:
        yield A
    else:
        for i in range(1,n+1):
            for hp in _heap_perm_(n-1, A):
                yield hp
            if (n % 2) == 1:
                j = 1
            else:
                j = i
            swap = A[j-1]
            A[j-1] = A[n-1]
            A[n-1] = swap

This produces permutations in the following order (for input of [0,1,2,3]):

0,1,2,3
1,0,2,3
2,0,1,3
0,2,1,3
1,2,0,3
2,1,0,3
3,1,2,0

and so on

This seems fine until that last step, which isn't a swap of one pair.

What am I doing wrong?

Answer 1

Historical prologue

The Wikipedia article on Heap's algorithm has been corrected since this answer was written but you can see the version referred to by the question and answer in Wikipedia's change history.

There's nothing wrong with your code (algorithmically), if you intended to implement the Wikipedia pseudocode. You have successfully implemented the algorithm presented.

However, the algorithm presented is not Heap's algorithm, and it does not guarantee that successive permutations will be the result of a single interchange. As you can see in the Wikipedia page, there are times when multiple swaps occur between generated permutations. See for example the lines

CBAD
DBCA

which have three swaps between them (one of the swaps is a no-op).

This is precisely the output from your code, which is not surprising since your code is an accurate implementation of the algorithm presented.

The correct implementation of Heap's algorithm, and the source of the error

Interestingly, the pseudocode is basically derived from Sedgewick's talk slides (reference 3 on the Wikipedia page), which does not correspond to the list of permutations on the immediately preceding slide. I did some sleuthing around, and found many copies of this incorrect pseudo-code, enough to start to doubt my analysis.

Fortunately, I also looked at Heap's short paper (reference 1 on the Wikipedia page), which is reasonably clear. What he says is: (emphasis added)

The program uses the same general method … i.e. for n objects, first permute the first (n—1) objects and then interchange the contents of the nth cell with those of a specified cell. In this method this specified cell is always the first cell if n is odd, but if n is even, it is numbered consecutively from 1 to (n−1).

The problem is that the recursive function as presented always does a swap before returning (unless N is 1). So it actually does swaps consecutively from 1 to n, not (n−1) as Heap specifies. Consequently, when (for example) the function is called with N==3, there will be two swaps at the end of the call before the next yield: one from the end of the N==2 call, and the other one from the loop with i==N. If if N==4, there will be three swaps, and so on. (Some of these will be no-ops, though.)

The last swap is incorrect: The algorithm should do swaps between recursive calls, not after them; it should never swap with i==N.

So this should work:

def _heap_perm_(n, A):
    if n == 1: yield A
    else:
        for i in range(n-1):
            for hp in _heap_perm_(n-1, A): yield hp
            j = 0 if (n % 2) == 1 else i
            A[j],A[n-1] = A[n-1],A[j]
        for hp in _heap_perm_(n-1, A): yield hp

Sedgewick's original paper

I found a copy of Sedgewick's 1977 paper, (the link given by wikipedia is paywalled, sadly), and was delighted to find that the algorithm he presents there is correct. However, it uses a looping control structure (credited to Donald Knuth) which might seem foreign to Python or C programmers: a mid-loop test:

Algorithm 1:

  procedure permutations (N);
      begin c:=1;
          loop:
              if N>2 then permutations(N-1)
              endif;
          while c<N:
              # Sedgewick uses :=: as the swap operator
              # Also see note below
              if N odd then P[N]:=:P[1] else P[N]:=:P[c] endif;
              c:=c+l
          repeat
       end;

Note: The swap line is taken from page 141, where Sedgewick explains how to modify the original version of Algorithm 1 (on page 140) to match Heap's algorithm. Aside from that line, the rest of the Algorithm is unchanged. Several variants are presented.