A recent post asks how to create a Maple permutation iterator, that is, a procedure that, on successive calls, iterates through each permutation of a given input. I suggested a routine that solved the problem, however, I wasn't satisfied with it. It was slower than it should be. Later I suggested an improvement. Here is another improvement. It uses the same algorithm (algorithm L, from The Art of Computer Programming, vol. 4, fasicle 2, by Donald E. Knuth) if the input has repeated elements, but uses a different method, algorithm T, ibid., if the input consists of distinct elements. The sequences of permutations for the two algorithms differ, the second algorithm uses the initial order for the first element and consecutive outputs differ by one transposition.
The output is returned as a one-dimensional rtable. This makes computing all permutations of reasonable sized lists efficient. Note that iterating through all permutations of, say, a 10 element list of distinct items is substantially faster and, more significantly, less memory hungry, than attempting to instantiate all elements using combinat[permute].
Permute := module()
export ModuleApply, Permute, Compile;
option load=Compile;
local AlgL, AlgT, compiled;
Compile := proc()
if compiled=true then
error "procedures are already compiled";
end if;
AlgL := Compiler:-Compile(AlgL);
AlgT := Compiler:-Compile(AlgT);
compiled := true;
NULL;
end proc;
ModuleApply := Permute;
Permute := proc(L :: Or(set,list,posint,range(integer)))
local u,n,DataType;
if L :: 'posint' then return procname([seq(u, u = 1..L)]);
elif L :: 'range' then return procname([seq(u, u = L)]);
end if;
DataType := 'datatype'='integer'[iquo(kernelopts('wordsize'),8)];
n := nops(L);
if n > 1 and (L :: set or n = nops(convert(L,'set'))) then
# elements are distinct
module()
export nextvalue, finished;
local complete,cnt,T,V,N;
V := rtable(convert(L,'list'));
N := n!;
T := rtable(1..N-1, DataType);
AlgT(T,n);
cnt := 0;
complete := evalb(n=0);
finished := () -> complete;
nextvalue := proc()
local j,k;
if complete then
error "attempt to call nextvalue on finished iterator";
end if;
if cnt<>0 then
j := T[cnt];
k := j+1;
(V[j], V[k]) := (V[k],V[j]);
end if;
cnt := cnt+1;
if cnt = N then
complete := true;
end if;
return V;
end proc;
end module;
else
# repeated elements
module()
export nextvalue, finished;
local complete,cnt,P,A,V;
A := rtable(1..n, sort(map(addressof,convert(L,'list')))
, DataType);
V := rtable(1..n, i -> pointto(A[i]));
cnt := 0;
P := rtable(1..n,1..2, DataType);
complete := evalb(n=0);
finished := () -> complete;
nextvalue := proc()
local i,j,k;
if complete then
error "attempt to call nextvalue on finished iterator";
end if;
# permute V per previous computation
for i to cnt do
(j,k) := (P[i,1],P[i,2]);
(V[j], V[k]) := (V[k],V[j]);
end do;
cnt := AlgL(A,P,n);
if cnt = 0 then complete := true end if;
return V;
end proc;
end module;
end if
end proc:
# This uses Algorithm L, pp 39-40, from
# "The Art of Computer Programming," vol 4, fasicle 2,
# "Generating all Tuples and Permutations", Donald Knuth
AlgL := proc(A :: Array(datatype=integer[BytesPerWord])
, P :: Array(datatype = integer[BytesPerWord])
, n :: posint
)
local j,k,l,cnt;
for j from n-1 to 1 by -1 do
if A[j+1] > A[j] then break; end if;
end do;
if j=0 then
return 0;
end if;
for l from n by -1 do
if A[l] > A[j] then break; end if;
end do;
(A[j],A[l]) := (A[l],A[j]);
cnt := 1;
P[1,1] := j;
P[1,2] := l;
l := n;
for k from j+1 do
if k >=l then break; end if;
(A[k],A[l]) := (A[l],A[k]);
cnt := cnt+1;
P[cnt,1] := k;
P[cnt,2] := l;
l := l-1;
end do;
return cnt;
end proc;
AlgT := proc(t :: Array(datatype=integer[BytesPerWord])
, n :: posint
)
local N,d,j,k,m;
N := n!;
d := N/2;
t[d] := 1;
k := N;
m := 2;
do
if k >= N then
if m=n then break; end if;
m := m+1;
d := d/m;
k := 0;
end if;
k := k+d;
j := m-1;
while j > 0 do
t[k] := j;
j := j-1;
k := k+d;
end do;
t[k] := t[k]+1;
k := k+d;
while j < m-1 do
j := j+1;
t[k] := j;
k := k+d;
end do;
end do;
return NULL;
end proc;
use bytesperword = iquo(kernelopts('wordsize'),8) in
AlgL := subs('BytesPerWord'=bytesperword, eval(AlgL));
AlgT := subs('BytesPerWord'=bytesperword, eval(AlgT));
end use;
end module:
LibraryTools:-Save(Permute,"Permute.mla");
Here is how it can be used.
libname := libname,currentdir():
#read "Permute.mpl";
# if you read it (rather than load from an archive),
# then you need to call the Compile export if you want to use the compiled procedures.
# Permute:-Compile();
N := 10:
P := Permute(N):
time(proc() while not P:-finished() do P:-nextvalue() end do end proc());
66.916
Comments
Combstruct improvement
Fantastic!
It looks as if all combstruct could be improved this way?
Alec
PS I remember times (70s, in Russia), when writing 2 lines of code per day (average) was considered being a very good performance :)
BytesPerWord
Hi Joe, is BytesPerWord the same as kernelopts(wordsize)/8 ?
You might prefer to make that integer datatype's width be dynamic in your code above, rather than hardcoded at preprocessor (ie. read) time with a $define.
acer
okay
Yes. I'll modify it. Doing so is slightly tricky and might be instructive. Part of the fun is doing this in a way that the resulting module can be saved to a library that is usable with either 32 or 64 bit Maple (that is, the byte-size is determined at run time and not load time). I think I got right. Note that if you read the file rather than loading it from a library you'll have to call the Compile export. Normally one could call the load procedure just before the end of the module assignment, however, that doesn't work here because when saving module the procedures cannot be compiled.
Actually I see I didn't get it quite right. The substitution in the final use statement at the end of the module assignment occurs when the module is assigned, so the type is set when the library is created. I originally had this in the procedure that is called when the module is loaded (here Compile), however, I was trying to modify it so that it would work for those users not accustomed to creating and loading libraries. Is there a nice way to handle that? I suppose I could put code into Permute to force the reassignment (and compilation), if not done, but that seems ugly.
Faster than next_permutation
I just compared it with next_permutation from STL, and Permute is more than 3 times faster!
Here are the details. First, I made a permute.dll from the cpp file with the following code,
#include <algorithm>
using namespace std;
__declspec(dllexport) BOOL np(int *a, int n){
return next_permutation(a,a+n);}
Then, in Maple,
First check that it is working,
A:=Array([1,2,3],datatype=integer[4]); A := [1, 2, 3] next_permutation(A,3); true A; [1, 3, 2]Now, time comparison
N := 10: P := Permute(N): time(proc() while not P:-finished() do P:-nextvalue() end do end proc()); 27.346 A:=Array([$1..10],datatype=integer[4]); A := [1, 2, 3, 4, 5, 6, 7, 8, 9, 10] time(proc() while next_permutation(A,10) do od end()); 90.801As I said already, fantastic!
Alec