<

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

Please Wait...