Heap’s algorithm permutation generator

Historical prologue

The Wikipedia article on Heap’s algorithm has been corrected, defaced and corrected again several times since this answer was written. The version referred to by the question and original answer was incorrect; you can see it in Wikipedia’s change history. The current version may or may not be correct; at the time of this edit (March 2022), the page contained both correct and incorrect versions.

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 the slides from a talk Sedgewick gave in 2002 (reference 3 on the Wikipedia page, also currently available on Sedgewick’s home page). The incorrect pseudocode is on slide 13, and it is clearly wrong since it does not correspond to the useful graphic showing the working of the algorithm on the immediately preceding slide. I did some sleuthing around, and found many copies of this incorrect pseudocode, 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

Note: The above function was written to mimic the function of the same name in the original question, which performs the permutation in-place. Doing the permutation in place saves the cost of copying every permutation returned, making the full generation O(n!) (or O(1) per permutation generated) instead of O(n·n!) (or O(n) per permutation generated). That’s fine if you’re going to process each permutation immediately, but confusing if your plan is to keep them around, since the next call to the generator will modify the previously generated list. In that case, you might want to call the function as

for perm in map(tuple, heap_perm(n, A)):
    # Now each perm is independent of the previous one

Or collect them all into a gigantic list (NOT recommended!!; the lists are huge) with something like perms = list(map(tuple, heap_perm(n, A))).

Sedgewick’s original paper

When I found Sedgewick’s 1977 paper [see note 1], I 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.

---

Notes:

1. When I wrote this answer, the only public link to the paper was a paywalled URL in a reference on the Wikipedia page. It is now available on Sedgewick’s home page.