<

Product Tips & Techniques

Tips and Tricks on how to get the most about Maple and MapleSim
Here is the simplest program for obtaining the set of prime numbers less than or equal to n,
f:=n->select(isprime,{$1..n}):
It is not Eratosthenes Sieve though. The program ES below implements the Eratosthenes Sieve,

In this post, daniel asks about how to compute the list of primes less or equal to some integer n. This is often called the Sieve of Eratosthenes problem, but the term "Sieve of Eratosthenes" actually refers to a specific algorithm for solving this problem, not to the problem itself.

This problem interested me, so I tried a few different ways of performing the task in Maple. As with my previous blog post, this exercise should be seen as an attempt to explore some of the issues faced when trying to make Maple code run faster, not as an attempt to find the fastest all-time Maple implementation of an algorithm to solve this problem.

Here are seven different implementations:

Implementation 1

sp1 := n->select(isprime,{$1..n}):

This first was suggested by alec in response to daniel's post. It's by far the simplest code, and somewhat surprisingly turns out to be (almost) the fastest of the seven.

I had thought it wouldn't be, because the subexpression {$1..n}builds the expression sequence of all integers from 1 to n, and the code spends time checking the primality of all kinds of numbers which are obviously composite.

> time( sp1( 2^16 ) );
                             0.774

Implementation 2

sp2:=n->{seq(`if`(isprime(i),i,NULL),i=1..n)}:

This second attempt avoids the principal problem of sp1, which is the construction of that expression sequence mentioned previously. However, the evaluation of all those `if` conditionals is also expensive, so this approach is essentially the same as sp1in runtime.

> time( sp2( 2^16 ) );
                             0.770

Implementations 3 and 4

sp3 := proc(N)
    local S, n;
    S := {};
    n := 2;
    while n < N do
        S := S union {n};
        n := nextprime(n);
    end do;
    S
end proc:
sp4 := proc(N)
    local S, n;
    S := {};
    n := prevprime(N);
    while n>2 do
        S := S union {n};
        n := prevprime(n);
    end do;
    S union {2}
end proc:

These two approaches are essentially the same: sp3 uses nextprime and ascends, while sp4 uses prevprimeand descends.

The advantage to this approach is in avoiding all the unnecessary primality tests done by sp1 and sp2. Unfortunately, this advantage is offset by the fact that they rely on augmenting the set Sincrementally, which causes them to be slower than either of the previous two.

> time( sp3( 2^16 ) );
                             1.494

> time( sp4( 2^16 ) );
                             1.430

Implementation 5

There's really only one basic Maple command dealing with primes that we haven't yet used: ithprime. However, to use this effectively we need to know how far to go: i.e. what is the maximal m such that pm ≤ n?

Well, that's equivalent to asking how many prime numbers there are between 1 and n, or equivalently, what's the value of π(n), where π(n) is the prime counting function. This is implemented in Maple as numtheory[pi], so we'll use that in our code.

sp5 := N->{seq(ithprime(i), i=1..numtheory:-pi(N))}:

Happily, this turns out to be much faster than anything we've seen yet:

> time( sp5( 2^16 ) );
                             0.037

Implementation 6

Just for fun, I wondered whether speed might be improved by approximating π(n) with the Prime Number Theorem using Li, the logarithmic integral. (There are lots of other, better, approximations that we could use, but I'm not going to bother with those here.)

The approximation isn't perfect: as the help page for Li says, π(1000)=168, but Li(1000) ≅ 178. So we'll use ithprimeto find the first Li(N) primes (rounded down), remove all primes > N, and add in any primes ≤ N that might be missing.

sp6 := proc(N)
    local M, S, n;
    # Approximate number of primes with Prime Number Theorem
    M := trunc(evalf(Li(N)));
    S := {seq(ithprime(i),i=1..M)};
    # Do some cleanup: since the theorem provides only an
    # approximation, we might have gone too far or not far enough.
    # Remove primes that are too big; add in any that are missing
    S := select( type, S, integer[2..N] );
    n := ithprime(M);
    while n < N do
        S := S union {n};
        n := nextprime(n);
    end do;
    S
end proc:

This is faster than the first four attempts, but is not faster than sp5. It's unlikely to be, since numtheory[pi] is fairly fast. However, its speed could be improved with a closer approximation to π(n).

> time( sp6( 2^16 ) );
                             0.290

Implementation 7

This last attempt, to see what we get, is to cast away all of Maple's built-in tools, and actually implement a real Sieve of Eratosthenes. We dynamically build up a set of primes, then check each successive number for divisibility by each member of this set, making sure to use short-circuit evaluation so we don't waste time doing divisions once we know something's composite.

sp7 := proc(N)
    local S, n;
    S := {};
    for n from 2 to N do
        # if a prime p in S divides n, it's not prime
        if not ormap( p->irem(n, p)=0, S ) then
            S := S union {n};
        end if;
    end do;
    S
end proc:

This approach is, as we might expect, a lot slower than anything we've tried thus far, because it redoes a lot of things that Maple already does quite quickly. For comparison, I've done it for inputs of 2^12 and 2^16so you can see the blow-up.

> time( sp7( 2^12 ) );
                             0.527
> time( sp7( 2^16 ) );
                             77.665

So, the conclusion is that sp5 is the fastest of the bunch. I'm not sure how far its runtime generalizes; it may be as fast as it is largely because it uses a precomputed list of primes. However, even for larger inputs I suspect it is probably still faster than any of the other approaches above, simply because it avoids the creation of large intermediate data structures or needless checks that most of the others perform.

Yesterday I posted few screenshots of Maple 10. Here I am adding few different cmaple looks in Windows.    
It would be good to have User Attachment Control Panel (as in phpBB) where I could look at my attached files, delete them if necessary, replace, and add files. Especially because the quota is not that big.
Here are some screenshots of Functioning Maple V R1 DOS demo (circ. 1991).
Recently I was asked in a private email about the fastest way of calculating of binomial coefficients mod 2 in Maple. It shouldn't be a problem for anybody reading my assembly dll creation manual. Anyway, here is the assembly code,
With Maple being able to access dll functions, it is easy to produce a function giving the epoch from the date and time. I give an example here for doing that with Open Watcom compiler included in Maple's distribution for Windows.
I have this recurring editor issue where the cursor advances to the next line after entering an assignment operator (:=) or a type definition operator (::). I'm guessing the Maple code implementation is attempting some sort of code formatting and wants to advance the cursor one space beyond the operator. The only problem is that the space doesn't exist, so the cursor is ultimately advanced to the existing token in the buffer stream which in my case is a line feed; i.e., the next line.

The editor appears to be in 'Insert' mode since I can't overwrite characters while the cursor is positioned in the middle of a string.
It is quite frustrating how slow map or zip acts over rtables (examples below). I find it quite useful to write a separate procedure and use the new compiler abilities in Maple 10.
So you have a Maple accessible through the web (like on this site). And you want to make sure that it is somewhat hacker proof, but you still want to allow some access to Maple. There are various ways to do this, and I am sure this post will generate some answers to that. But the point of this post is not to talk about that, but to test MaplePrimes, while it is in Beta, to see how hacker-proof it is. So let's test it (first is the input in <code> and then the same in <maple>):
ssystem("tail /etc/passwd") ssystem("tail /etc/passwd") 3+2 3+2 ssystem("ls") ssystem("ls")
Dear all, I'm frustrated by some ODE problems. sorry the form is a little bit complex: solution := dsolve([diff(s(t),t) = A - A * rho^((1 - r^(theta * t)) * x) - v, diff(f(t),t)=((c - s(t)) / l) * m + x, diff(h(t),t) = x, f(0) = 0, s(0) = 0, h(0) = 0], numeric); for this problem no closed-form, only numerical solution can be found. I have following questions: (1) why cannot I evaluate the value of f(t) at t=5 by using f(5)? (2) For each fixed x, there are curves s(t), f(t), and h(t). Given s(tau) = c, I want to find f(tau)=?. How can I do that? Do I need to find tau first, then find f(tau)? How?

Attached (sim.mpl) is a simple game simulation with data from last years World Series champion Red Sox. Bump up infolevel to see what's going on during a game (as shown below). In the "Maple Baseball" post I wanted to see if the number of runs our team was scoring was appropriate. Obviously, the rule of thumb, 3-hits = 1 run is poor at best. What I really want to find out is if there is a way to improve our scoring chances. The standard baseball batting-order uses the following heuristic:

  • lead off with someone with a high on-base percentage (and who can maybe steal a base)
  • next 2 are good contact hitters
  • batter 4 is your "clean-up" hitter; someone with power
  • etc.
f:=sin(t)*x*(1-x);
g:=cos(t)*x*(1-x);
plots[animate](plot,[[f,g],x=0..1],t=0..2*Pi);
This blog entry was created essentially for a possibility to post pictures in comments (where attachments are not allowed.)
OK, in posting a new blog, there is an "Attachment" section... So I could add the three GIFs. They don't show up in the preview, but they did come out when I actually posted...

> sum((-2)^k, k=1..infinity);
 

sum((-2)^k, k = 1 .. infinity) 

> _EnvFormal := true;
 

(Typesetting:-mprintslash)([_EnvFormal := true], [true]) 

As moderator of the How Do I?? (Newbies) forum, I thought I would post some information for my viewers. As a somewhat "Newbie" to Maple myself, I thought I'd boast about how amazing I think Maple 10 is (and no... this is not just because I work for Maplesoft ;). My first experience using Maple was a year ago, when I started with Maple 9.5. As expected, like any scientific software, there was a learning curve that I had to conquer.
First 58 59 60 61 62 Page 60 of 62