Faster Permutations

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

alec's picture

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 :)

acer's picture

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.

alec's picture

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,

next_permutation:=define_external(
'`?np@@YGHPAHH@Z`',
'a'::ARRAY(datatype=integer[4]),
'n'::integer[4],
'RETURN'::boolean[4],
'LIB'="permute.dll"):

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.801

As I said already, fantastic!

Alec

Comment viewing options

Select your preferred way to display the comments and click "Save settings" to activate your changes.
}