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?
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.
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))).
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.