Generate Permutations using Heap's Algorithm

Question

I am trying to generate all permutations for an array using the Heap's algorithm that I found in wikipedia.

This is what I tried so far:

n <- 3
A <- c(1, 2, 3)
perm <- function(n, A) {
  if (n == 1)
    print(perm)
  for (i in length(A))
    perm(n, A - 1) 
  if (A %% 2 == 1) {
    swap(A[i], A[n - 1])
  } else {
    swap(A[0], A[n - 1])
  } 
}
perm(3, A)

The output is not showing, and it would be great to receive some help.

Answer 1

9In addition to some misuse of names, there are four fundamental issues with the provided code:

1. The implementation of Heap's algorithm in Wikipedia assumes 0-based indexing, while R's vectors are 1-based. So A[0] needs to be translated to A[1] (that is, the first element in the array), while A[n-1] needs to be translated to A[n] (which is the last element in the n-prefix of A).

2. As is surprisingly common, the code does not follow the iterative structure of the Wikipedia code. Stripped to its essence, the correct code is (changed to 1-based indexing):

   define heap(A, n):
     if n == 1: process A
     else:
       for i from 1 to n - 1:    ## Note the n - 1
         heap(A, n - 1)          ## Recurse
         swap A[n] with ...      ## Swap
       end for
       heap(A, n - 1)            ## Recurse
     end if
   

   The important aspect of this loop is that the swap is done in between recursions, neither before nor after. So there are n recursive calls (assuming n > 1) but only n - 1 swaps. The consequence is that there is exactly one swap between two successive permutations, which is the essence of the "Gray code" requirement which Heap's algorithm satisfies.

3. The various pseudocode implementations (and the concrete implementations in Python) assume that the array argument is actually passed by reference, so that swaps in the recursive calls are visible to the caller. But R passes arrays by value, so every recursive call has its own instance of the array, and the modifications are not reflected in the caller's array.

4. The computation of the swap is as follows:

   • If n is even, A[n] is swapped with A[1].
   • If n is odd, A[n] is swapped with A[i].

   (Remember that the loop limit is n - 1, not n, and n is always greater than 1. So A[n] is never swapped with itself.)

   The provided code gets this backwards.

To get around the issue with making a copy of the array at every recursive call, we can use a closure containing the array; the only argument for the recursive call is then the prefix length. But note that modifying the array in the closure's environment requires that we use the <<- operator.

heap_perm <- function(A) {
  h <- function(n) {
    if (n == 1) {
      print(A)
    } else {
      h(n - 1)
      for (i in 1:(n - 1)) {
        if (n %% 2 == 1) pivot <- 1 else pivot <- i
        tmp <- A[n]; A[n] <<- A[pivot]; A[pivot] <<- tmp
        h(n - 1)
      }
    }
  }
  
  h(length(A))
}

When testing implementations of Heap's algorithm, it's important to use a vector of at least 4 elements. (A good test would use a larger vector, but four is a good start.) When examining the output, you need to check that:

• The correct number of permutations are produced (24 in the case of a vector of 4 elements);
• Each generated permutation is unique;
• Two consecutive permutations differ only in a single swap.

Verifying all three of these conditions will avoid the most common bugs in implementations of Heap's algorithm.

Here's the test output of the above R program:

> heap_perm(0:3)
[1] 0 1 2 3
[1] 1 0 2 3
[1] 2 0 1 3
[1] 0 2 1 3
[1] 1 2 0 3
[1] 2 1 0 3
[1] 3 1 0 2
[1] 1 3 0 2
[1] 0 3 1 2
[1] 3 0 1 2
[1] 1 0 3 2
[1] 0 1 3 2
[1] 0 2 3 1
[1] 2 0 3 1
[1] 3 0 2 1
[1] 0 3 2 1
[1] 2 3 0 1
[1] 3 2 0 1
[1] 3 2 1 0
[1] 2 3 1 0
[1] 1 3 2 0
[1] 3 1 2 0
[1] 2 1 3 0
[1] 1 2 3 0

Answer 2

To implement Heap's algorithm in R, here is a solution that outputs the result in a matrix for further processing:

h <- function(s, n, e, ...) {
  ## processes Heap's algorithm
  if (n == 1L) {
    e$i <- e$i + 1L
    e$r[e$i, ] <- e$s
  } else {
    h(s, n - 1L, e)
    for (i in 1:(n - 1L)) {
      if (n %% 2L == 0L) {
        j <- i
      } else {
        j <- 1L
      }
      tmp <- e$s[n]
      e$s[n] <- e$s[j]
      e$s[j] <- tmp
      h(s, n - 1L, e)
    } 
  }
}

hperm <- function(s) {
  ## wrapper
  e <- environment()
  l <- length(s)
  i <- 0L
  e$r <- matrix(NA, factorial(l), l)
  h(s, l, e, i)
  return(e$r)
}
    

Gives:

A <- 0:3

hperm(A)
#       [,1] [,2] [,3] [,4]
#  [1,]    0    1    2    3
#  [2,]    1    0    2    3
#  [3,]    2    0    1    3
#  [4,]    0    2    1    3
#  [5,]    1    2    0    3
#  [6,]    2    1    0    3
#  [7,]    3    1    0    2
#  [8,]    1    3    0    2
#  [9,]    0    3    1    2
# [10,]    3    0    1    2
# [11,]    1    0    3    2
# [12,]    0    1    3    2
# [13,]    0    2    3    1
# [14,]    2    0    3    1
# [15,]    3    0    2    1
# [16,]    0    3    2    1
# [17,]    2    3    0    1
# [18,]    3    2    0    1
# [19,]    3    2    1    0
# [20,]    2    3    1    0
# [21,]    1    3    2    0
# [22,]    3    1    2    0
# [23,]    2    1    3    0
# [24,]    1    2    3    0

See also: This Python code of Heap's algorithm on Stack Overflow.