Joe Riel

9530 Reputation

23 Badges

20 years, 27 days

MaplePrimes Activity


These are answers submitted by Joe Riel

N.B. This has been completely rewritten.

On rethinking this, I realized that this can be improved by returning a Vector rather than a list.  That has two advantages: it saves memory (each new output is not inserted into Maple's Simplification table) and it is faster because we don't have to process every element (with a call to pointto).  It also permits using the Maple compiler to handle the core of algL.  Here is the rewrite. 

Note that I used a preprocessor directive, BytesPerWord, so that the source code could conceivably be used with 64 bit Maple; however, the use of a preprocessor directive means it will only work if run through command-line maple (my usual mode of working).  You can always manually replace the few occurrences of BytesPerWord with 4 for 32 bit Maple or 8 for 64 bit Maple.  I couldn't think of a nice way to declare the Arrays in ComputePerm so that using the preprocessor directive was not required.

$define BytesPerWord 4

Permute := module()
export ModuleApply, Permute;
local ComputePerm;

    ModuleApply := Permute;
    
    Permute := proc(L :: Or(set,list,posint,range(integer)))
    local u,n,A,V;
        
        if L :: 'posint' then return procname([seq(u, u = 1..L)]);
        elif L :: 'range' then return procname([seq(u, u = L)]);
        end if;
        
        n := nops(L);
        A := rtable(1..n, sort(map(addressof,convert(L,'list'))), 'datatype=integer[BytesPerWord]');
        V := rtable(1..n, i -> pointto(A[i]));
        
        module()
        export nextvalue, finished;
        local complete,cnt,P;
            
            cnt := 0;
            P := rtable(1..n,1..2, 'datatype=integer[BytesPerWord]');
            
            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 := ComputePerm(A,P,n);
                if cnt = 0 then complete := true end if;
                return V;
            end proc;
        end module;
    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
                
    ComputePerm := 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;
    
    ComputePerm := Compiler:-Compile(ComputePerm);
    
end module:

 

P := Permute([1,2,3]):
while not P:-finished() do printf("%d\n",P:-nextvalue()); end do;
1 2 3
1 3 2
2 1 3
2 3 1
3 1 2
3 2 1


N := 10:
P := Permute([seq(1..N)]):
time(proc() while not P:-finished() do P:-nextvalue() end do end proc());

                                                     135.0

 

We resolved this offline. In the unlikely case this is of interest to any one else, the problem was that 'nil', indicating that no maple initialization file is specified, was missing from maplev-executable-alist (the missing 'nil' exists in the assignment shown in the original post).  Possibly an issue using customize-group, but I wasn't able to reproduce it. 

I'm curious, does it work on Maple 9 with a range or integer input?  I suspect it may not, I think that seq(1..3) wasn't added until Maple10, but am not sure and don't have Maple 9 quickly accessible. It would be easy enough for you to fix that, just modify the start, where ranges are converted to lists of integers (if you care).

Alvaro,

Do you use GNU Emacs or Xemacs?  I'll assume the former (which is almost required).  Do other parts of the mode appear to work (highlighting, auto-indentation, etc)?   Try running the emacs command maplev-mint-buffer directly:

M-x maplev-mint-buffer

where M-x means use the meta (alt) modifier key while pressing x.  If that works, then check what key it is mapped to (in a maplev buffer):

C-h w maplev-mint-buffer

 

Note that from Maple's perspective, it makes little sense to ask for permutations of a set.  All permutations of an unordered set are equal.  With that in mind, the following should do what you want.  It accepts either a set, a list, a posint, or a range of integers as input, but always returns lists.  Being somewhat lazy this morning, I originally wrote it to handle only integers, then to generalize it used the hackware procedures addressof and pointto.  It returns a module with two exports.  In that regard it is slightly different from combinat[cartprod], that is, the export finished is a procedure not a boolean value:

GeneratePermutations := proc(L :: Or(sequential,posint,range(integer)))

    if L :: 'posint' then
        return procname([seq(1..L)]);
    elif L :: 'range(integer)' then
        return procname([seq(L)]);
    end if;

    module()
    export nextvalue, finished;
    local complete,n,A;

        n := nops(L);
        A := rtable(1..n,sort(map(addressof,convert(L,'list'))));

        complete := `if`(L=[], true, false);
        finished := () -> complete;

        # This uses Algorithm L, pp 39-40, from
        # "The Art of Computer Programming," vol 4, fasicle 2,
        # "Generating all Tuples and Permutations", Donald Knuth

        nextvalue := proc()
        local a,j,k,l;

            if complete then
                error "attempt to call nextvalue on finished iterator";
            end if;

            a := map(pointto,convert(A,'list'));

            # update A

            for j from n-1 to 1 by -1 while A[j+1] <= A[j] do end do;
            if j=0 then
                complete := true;
                return a;
            end if;

            for l from n by -1 while A[l] <= A[j] do end do;
            (A[j],A[l]) := (A[l],A[j]);

            # reverse A[j=1]..A[n]
            l := n;
            for k from j+1 while k < l do
                (A[k],A[l]) := (A[l],A[k]);
                l := l-1;
            end do;

            return a;
        end proc;
    end module;
end proc:

Here is an example

P := GeneratePermutations([a,[1,2],[1,2]]):
while not P:-finished() do P:-nextvalue(); end do;
                              [a, [1, 2], [1, 2]]
                              [[1, 2], a, [1, 2]]
                              [[1, 2], [1, 2], a]

P := GeneratePermutations(2..4):
while not P:-finished() do P:-nextvalue(); end do;
                                   [2, 3, 4]
                                   [2, 4, 3]
                                   [3, 2, 4]
                                   [3, 4, 2]
                                   [4, 2, 3]
                                   [4, 3, 2]

You might want to look at the CodeTools[Profiling] subpackage.  Look at the help pages for Profile and PrintProfiles. Profiling is useful for  determining the sections that deserves the most attention.

See the help page for LinearAlgebra[MatrixInverse].  It can compute the Penrose pseudo-inverse.  The LinearAlgebra package also has IndentityMatrix and ZeroMatrix procedures.  See LinearAlgebra,Details for a table of all procedures.

You are comparing solutions to different problems.  Note that the original post had a typo (read down a few replies to see the correction).  His original description had the exponent as 3/2.  It should have been 2/3.  My original response (with the A and B) was for the 3/2 exponent.  My later response was for the 2/3 exponent.

You don't need to put it in that form.  Just do

numer(x+y/z);
                        x*z + y

Well, you could always cheat and use identify:

len := VectorCalculus:-ArcLength(<x,(x/2)^(2/3)>, x=0..2):
identify(evalf[50](len));                                 
                                                    1/2
                                               20 10
                                               -------- - 2/27
                                                  27

This works,

len := VectorCalculus:-ArcLength(<x,(x/2)^(2/3)>, x=0..2,inert):
IntegrationTools:-Change(len, y=x^(2/3)):
value(%);
                                                    1/2
                                               20 10
                                               -------- - 2/27
                                                  27

#1:  v is assigned, sequentially, each element off L.  v[2] is the second element of v.  For example,

L := [[a,f1],[b,f2],[c,f3]]: # (normally fk are numerical values, here I made them symbolic for demonstration.
tot := add(v[2], v in L);
                                 tot := f1 + f2 + f3

So, what it does is sum the frequencies in the list and assign the total to tot.  I'm not sure what you are implying with  your comment about just counting the total elements in the list.  The Shannon-Fano algorithm does not use  (half) the number of elements but rather (half) the sum of the frequencies as the criteria.

#2: Here I use k as an index. The reason for doing this rather than

for v in L while cumul < tot/2 do
    cumul := cumul + v[2];
end do

is that the purpose of this loop is to determine the index, k, so that the sum of the frequencies in elements 1 to k is approximately half the total.  Another way to do this, possibly clearer, is to use a separate counter rather than the loop counter:

split := 0:
cumul := 0:
for v in L while cumul < tot/2 do
    split := split+1;
    cumul := cumul + v[2];
end do;

A difference with that approach is that when the loop ends, split is one less than the k was in my original code. So you would have to adjust the ranges accordingly.

Note that in both implementations (my original and this trivial modification) the split is such that the sum of the frequencies of the elements in the first list is less than half the total.

If you want to see the code in action, call

debug(ShannonFano);
before making the call to ShannonFano.

Use the solve procedure.  Combine the two equations into a set, and pass them to solve.  For this problem, in which the coefficients are numerical, you can use fsolve, which is generally faster (not an issue here).  Note that in Maple multiplication is expressed with *.  For example:

eqs := {3*x + 4*y = 1,  y = 3}:
solve(eqs, {x,y});
fsolve(eqs, {x,y});

Something like this is closer to what you had in mind, I think.  I may have gotten the details of the implemenation wrong, but this should give you the idea.  Note that I return a nested list that encodes the tree structure.  Also note that I changed your compare to >, which requires less looping.

compare := proc(a::list,b::list)
    evalb(a[2] > b[2])
end proc:

ShannonFano := proc(L :: listlist)
local n,k,v,tot,cumul;
    n := nops(L);
    if n < 2 then
        return L[];
    end if;
    tot := add(v[2], v in L);
    cumul := 0;
    for k to nops(L) while cumul < tot/2 do
        cumul := cumul + L[k,2];
    end do;
    [procname(L[1..k-1]), procname(L[k..n])];
end proc:

listA :=[[a,41],[b,33],[c,23],[d,25],[e,100]]:
listB := sort(listA, compare);
            listB := [[e, 100], [a, 41], [b, 33], [d, 25], [c, 23]]

ShannonFano(listB);
             [[[e, 100], [a, 41]], [[[b, 33], [d, 25]], [c, 23]]]


Try

Re(exp(I*k*x)) assuming real;
                                   cos(k x)

A minor correction.  It is not true that explicit looping in Maple is necessarily inefficient.  In fact, for the particular task given, which consists of printing a column of primes, most of the time is likely to be spent in the print routine, if not in the primality testing. Note by the way (this is intended for the original poster) that printing, required here, is not something one ordinarily does with Maple.  It is usually more effective to return output (say as a list of primes) rather than to print it, since we can further manipulate the output but can do little but look at printed output.

To demonstrate my point about explicit looping not necessarily being inefficient, here are two procedures that return a list of primes, from m to n.  The second one, which is essentially what Doug described, is the superior procedure. It is short, clear, and does the job. However, it really is is only marginally faster than the first.  The call to time shows the elapsed duration required for each to return all primes from 1 to 100,000.

GetPrimes1 := proc(m,n)
local cnt,T,i;
    cnt := 0;
    T := table();
    for i from m to n do
        if isprime(i) then
            cnt := cnt+1;
            T[cnt] := i;
        end if;
    end do;
    convert(T,'list');
end proc:

GetPrimes2 := proc(m,n)
    select(isprime,[seq(m..n)]);
end proc:
N := 100\000:
time(GetPrimes1(1,N));
                                     1.880

time(GetPrimes2(1,N));
                                     1.824

 

First 93 94 95 96 97 98 99 Last Page 95 of 114