I'm really not very familiar with clojure. I think this is probably the longest clojure program I've written so far. But I guess it provides some sort of answer.
On the whole, functional enumerations of permutations are going to suffer from the O(n) cost of creating a new permutation vector on each iteration. Pretty well all of the common imperative solutions make some attempt to avoid this cost, and many of them can produce amortised O(1) complexity (for a single permutation). Heap's algorithm is interesting in an imperative environment largely because it guarantees to perform a single swap on each iteration, which is clearly a minimal total number of mutations. That's particularly useful if the goal is to perform some sort of aggregate computation over the permutation which can be incrementally computed from the previous value and the modifications.
Heap's algorithm is not the only algorithm which performs just a single swap to produce the next permutation. The possibly even more famous bell-ringers' algorithm (often called the Steiner-Johnson-Trotter algorithm) produces sequences in which consecutive permutations differ only by a swap of two adjacent elements. This could be even more valuable for updating of aggregate computations, but it is more complicated to figure out which two adjacent elements to swap at each iteration. (Indeed, although tables of plain changes for up to seven bells were produced several centuries ago, the precise algorithm used to produce these table has not, as far as I know been recorded. However, the results line up with the algorithm proposed independently in the 1960s by the three mathematicians after whom it is named.)
The plain changes algorithm was motivated in part by the desire to keep the various bell ringers' attention on the changes; if a bell stays in the same position in the change for too long, its ringer may become bored and lose their place in the sequence. (In fact, current change ringing sequences try even harder to avoid leaving a bell in the same position for too long, so they no longer use the STJ algorithm.)
For many combinatorial problems, though, the opposite criterion is desired. For example, it may be useful to generate the permutations in lexicographic order, which means that the first item retains its value for (n-1)! iterations. Lexicographic ordering can still be performed in amortised O(1), and it involves O(1) mutations at each iteration, but the number of mutations may be as great as n-1 (or n, if the permutation sequence is circular, since the last permutation in lexicographical order is the reverse of the first one).
If the vector to be permuted may have repeated elements, and only unique permutations are desired, then lexicographic order is far and away the easiest solution. Furthermore, it is very easy to describe the next-permutation algorithm for lexicographic ordering:
Start with the vector sorted from left to right.
At each iteration:
Find the shortest suffix which is not monotonically non-increasing. (In other words, find the last element for which some subsequent element is greater.) If the entire vector is monotonically non-increasing, then it is the last permutation and the process is done.
Sort that suffix by shifting its first element to the right, and then reverse it. (This can be done by a modified reversal algorithm with k/2 swaps where k is the length of the suffix, but doing it in two steps is conceptually simpler.)
I think that algorithm would be quite simple to implement in functional style.
But let's get back to Heap's algorithm. Heap developed his algorithm in 1963, about the same time as Trotter and Johnson independently published the change ringing algorithm. (Steinhaus had published his version of the algorithm in 1958, but since he was writing in Polish it was relatively unknown until it was translated into English in 1963.) Heap's paper was a simple description of his algorithm, without formal proof, and it probably would have remained in the dusty archives of computer science had it not been rediscovered by Robert Sedgewick a decade later. In 1977, Sedgewick wrote a long survey on permutation algorithms, in which he devoted quite a bit of attention to Heap's algorithm, including producing an optimised implementation in a kind of virtual machine code. He concluded that for typical hardware architectures (in 1977), "Heap's method will run faster than any other known method."
A couple of years later, he presented a public lecture on permutation algorithms, during which he repeated this claim, and Heap's algorithm became the permutation algorithm of choice for programmers with a performance fetish. Unfortunately, the presentation slides (which are rather more easily found online than any of the other sources for Heap's algorithm) had a minor error in the algorithm pseudocode (not present in the 1977 paper, which presents a number of variant implementations, all impeccable). And because the presentation slides are rather more accessible than any of the academic references, that particular error has plagued implementations of Heap's algorithm ever since. Ironically, the error causes the algorithm to both run more slowly and to generate an incorrect sequence in which consecutive permutations sometimes do not differ by a single swap.
I have a certain affinity to Heap's algorithm because a few years ago, a question appeared on StackOverflow asking for help with that algorithm; the poster had carefully implemented the algorithm based on pseudocode in Wikipedia, and discovered that their code did not work as expected. In the course of reviewing the code, I soon realised that the poster's code was an entirely accurate implementation of the Wikipedia write-up, and in the course of trying to validate Wikipedia, I tracked down the history of the error which I summarised in the preceding paragraph. A Wikipedian noticed the answer and fixed the Wikipedia page. (Constant vigilance is necessary, however; every couple of months some aspiring computer science expert compares the pseudocode on Wikipedia either with Sedgewick's buggy slide or with one of the many buggy implementations floating about the internet, and "corrects" the Wikipedia page.)
The beauty of Heap's algorithm is the simplicity of enumerating the swaps. The algorithm is similar to a lexicographic algorithm, in reverse, in that it always does a complete permutation cycle of a prefix of the elements before changing the next element. So the higher index of the swapped elements follows a kind of factorial variant on the ruler function, in which element i appears every i! iterations (except for the iterations which are multiples of a larger factorial). The lower index follows an only slightly more complex pattern: for prefixes with an odd number of elements, the lower index is always 0; for prefixes with an even number of elements, the lower index starts at 0 and is sequentially incremented. So, for example, in the case of a four element permutation sequence, the swaps are:
0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 2 0 0 0 0 0
1 2 1 2 1 3 1 2 1 2 1 3 1 2 1 2 1 3 1 2 1 2 1
(Note that the element with index 3 appears as the larger index every 3! == 6 iterations, and the element with index 2 appears as the larger index every 2! == 2 iterations, except the ones which have already been taken by index 3.)
This sequence of pairs of indices is very easy to produce recursively, but one really wants it as a sequence. There is a non-recursive algorithm which effectively involves maintaining an explicit stack, but since the published algorithms tend to mutate this stack, I chose to use the recursive algorithm. In many Scheme's, one would use some kind of continuation to turn the recursive algorithm into a sequence, but since clojure doesn't have continuations (as far as I could see), I did it with nested calls to mapcat, the same approach as is taken by clojure's tree-seq standard library function.
Once you have the list of pairs of indices to swap, they can be repetitively (and lazily) applied to permutations by using reductions, which generates the result of successively applying a sequence of values to a seed value, using an arbitrary binary function. You can see that in the last line, where reductions is applied to the sequence generated by swaps using the swap-v function to turn each successive permutation into the next one.
I hope that's enough narrative to explain the (possibly non-idiomatic) clojure program which I came up with:
(letfn [(even-swaps [n]
(if (= n 2) [[0 1]]
(drop 1 (mapcat (fn [i] (cons [(- i 1) (- n 1)]
(odd-swaps (- n 1))))
(range n)))))
(odd-swaps [n]
(drop 1 (mapcat (fn [i] (cons [0 (- n 1)]
(even-swaps (- n 1))))
(range n))))
(swaps [n]
(cond
(< n 2) '()
(even? n) (even-swaps n)
:else (odd-swaps n)))
(swap-v [v [i j]] (assoc v i (v j) j (v i)))]
(defn heap [sq]
(let [v (vec sq)]
(reductions swap-v v (swaps (count v))))))
Sample output:
user => (run! prn (heap '(a b c d)))
[a b c d]
[b a c d]
[c a b d]
[a c b d]
[b c a d]
[c b a d]
[d b a c]
[b d a c]
[a d b c]
[d a b c]
[b a d c]
[a b d c]
[a c d b]
[c a d b]
[d a c b]
[a d c b]
[c d a b]
[d c a b]
[d c b a]
[c d b a]
[b d c a]
[d b c a]
[c b d a]
[b c d a]
(defn swap-v [v i j] (assoc v i (v j) j (v i)))$\endgroup$