6
$\begingroup$

As Mathematica's implemented Permutations function is not compilable, I tried to write my very own Permutations implementation, called PermutationsNew, which I want to compile later on. Unfortunately the only implementation I came up with so far uses recursive programming and looks like this:

PermutationsInternal = Function[{list, listfixed}, If[Length[list] == 1,
                                Join[listfixed, list],
                                Table[PermutationsInternal[Drop[list, {i}],
                                      Append[listfixed, list[[i]]]], {i, Length[list]}] ]];
PermutationsNew[list_] := Flatten[PermutationsInternal[list, {}], Length[list] - 2]

The implementation seems to work fine, but because of the recursive structure, I cannot come up with any approach that can be compiled, because Mathematica has trouble compiling the lists encountered in PermutationsInternal.

Does anybody have an idea on how to compile the above approach or can offer me an idea for a different compilable solution?

Improve this question
$\endgroup$
6
  • $\begingroup$ What is your use case? With permutations there's a danger of running out of memory before running out of time ... $\endgroup$
    – Szabolcs
    Mar 29, 2014 at 17:17
  • $\begingroup$ I mean, Permutations returns a packed array, so it is memory efficient. ByteCount@Permutations[Range[11]] // AbsoluteTiming gives 1.8 seconds and 3.5 GB on my machine. There's clearly not enough memory to go to $n=12$ and 1.8 seconds is still not that slow for a one-time computation. $\endgroup$
    – Szabolcs
    Mar 29, 2014 at 17:22
  • $\begingroup$ memory usage should not be a problem, the lists I need to permute usually do not contain more than 8 items, but the permutation process would be part of a highly optimized and compiled function that in return has to be applied frequently in the computation process $\endgroup$
    – Wizard
    Mar 29, 2014 at 17:24
  • 1
    $\begingroup$ In that case how about precomputing permutations once, perm = Permutations@Range[4], then applying these precomputed permutations to your actual list, as in {a,b,c,d}[[#]]& /@ perm? This application can be done using a compiled function much more easily. If your list has the same integers as Range[1,n], then the permuted lists will be equivalent to perm, except in a different order. So this will probably make sense if your list has other numbers than {1,2,3,4} in some order. $\endgroup$
    – Szabolcs
    Mar 29, 2014 at 17:28
  • $\begingroup$ I think that's a good solution. I will try this out soon, but I do not see why it should not work. Thank you very much. $\endgroup$
    – Wizard
    Mar 29, 2014 at 20:07

2 Answers 2

Reset to default
9
$\begingroup$

As a proof of concept, I present here a straightforward implementation of Don Knuth's "Algorithm L" for generating permutations in lexicographic order:

cPermutations = Compile[{{list, _Integer, 1}},
                        Module[{n = Length[list], p = Sort[list], pbag, j, k},
                               pbag = Internal`Bag[Most[{1}]];
                               While[True,
                                     Internal`StuffBag[pbag, p, 1];
                                     j = n - 1;
                                     While[j > 0 && p[[j]] >= p[[j + 1]], j--];
                                     If[j == 0, Break[]];
                                     k = n;
                                     While[p[[j]] >= p[[k]], k--];
                                     p[[{j, k}]] = p[[{k, j}]];
                                     p = Take[p, j] ~Join~
                                         Take[p, {n, j + 1, -1}]];
                               Partition[Internal`BagPart[pbag, All], n]],
                        RuntimeOptions -> "Speed"];

It gives the same results as Permutations[]:

cPermutations[Range[9]] === Permutations[Range[9]]
   True

l = {1, 2, 2, 2, 3, 3, 5, 5, 7};
cPermutations[l] === Permutations[l]
   True

Unfortunately, this implementation is quite a lot slower than Permutations[]. There might be a few tricks to speed this up further, but considering that Permutations[] has been a built-in function since the first version, I am skeptical that it can be outdone by a compiled function like this one.

Improve this answer
$\endgroup$
4
$\begingroup$

This is not an answer, because unfortunately, recursion is not supported by the compiler. We can see this with a very simple example: the Fibonacci sequence. Recall that one can write a recursive pure function by using slot #0 to refer to the function itself. Hence:

With[{fibonacci = Function[Null, If[#1 == 1 || #1 == 2, 1, #0[#1 - 1] + #0[#1 - 2]]]},
 Compile[{{i, _Integer, 0}}, fibonacci[i]]
]

However, this input crashes the kernel. It takes some seconds for this to happen, so I expect it is caused by the attempts of the compiler to unroll the recursion, resulting in a stack overflow. Your function produces the same outcome if written in this way, but I felt that a more minimal example (without e.g. Append) would be helpful in demonstrating the source of the problem.

Improve this answer
$\endgroup$
2
  • $\begingroup$ "recursion is not supported by the compiler", this is not entirely true. If you execute this line twice, you get a recursive library call that avoids MainEvaluate: fib = Compile[{{n, _Integer}}, If[n == 1 || n == 2, 1, fib[n - 1] + fib[n - 2]], CompilationOptions -> {"InlineExternalDefinitions" -> True, "InlineCompiledFunctions" -> False}, CompilationTarget -> "C"] $\endgroup$
    – Marius Ladegård Meyer
    Aug 30, 2015 at 20:18
  • 1
    $\begingroup$ @MariusLadegårdMeyer it is true that one can use LibraryLink functions to get around the limitation, but I would argue that the existence of the opcode 42 mechanism does not equate to support for recursion per se. At that point it is easier to just write the function in C and not use the Mathematica compiler at all. $\endgroup$
    – Oleksandr R.
    Aug 30, 2015 at 22:49

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.