Items tagged with integer

Feed
Silvexil announced in his blog a Binary Arithmetic package. Here is my version of it.
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.

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

Chi^2 calculations above some "size" or "complexity", using Maple 9.5 and Global Optimization Toolbox (GOT), may produce after some time of calculation error messages like:

"Execution stopped: stack limit reached.
The kernel has been shut down. Further computation cannot be performed."

Seeking workarounds, I have looked for information at ?kernelopts for kernelopts(stacklimit), but it was not very useful:

"Limits may be raised or lowered. Maple limits may not be raised above any system defined hard limits. "

It would be convenient if the subscripted version of type/integer could handle infinity and -infinity. Then, to specify an integer greater than, say, 1, we could do type(i, integer[2..infinity]). Currently I handle this as type(i, And(integer,Range(1,infinity))) which is not as nice, particularly because it isn't clear that 1 is excluded. The drawback of doing this is that it implies that infinity is allowed. However, because infinity is not an integer, it seems reasonable that it would return fals
The IBM Research August 2005 Ponder This challenge is out. The attached 11 line Maple procedure solves it in just under 2.5 seconds. Don't look at my solution if you want to do this yourself.
Previously I described how to change the default zoom setting for the Maple gui by modifying the appropriate initialization file. Another useful setting to change is the default background color of the help browser. This is done by modifying, in the initialization file, the line HelpBGColor=. I set it to HelpBGColor=240 240 240, that gives a light gray background that is less harsh on my eyes. The three fields should be integers from 0 to 255; they correspond to the red, green, and blue components of the color.

A while ago, I was trying to determine the fastest way, using Maple, to compute a list of a particular sequence of integers defined by a recurrence relation. This isn't a new idea of course, but I think the exercise is a useful introduction into the various ways of coding algorithms in Maple, and why some may be better than others.

My original example was a bit esoteric, so for simplicity I've redone it using the standard example of a recurrence relation, the Fibonacci sequence. We'll fix N, the number of Fibonacci numbers to compute, at 40000.

> N := 40000:

Below, we'll try to find the fastest way, in Maple, to compute a list of these numbers from scratch. Note that the speed comparisons might be specific to this machine, OS, or Maple version, and even though the runtime of most of these implementations is about the same, the fastest method at N=40000 might not be the fastest one for larger N.

The first idea is the completely naive approach. We'll store the Fibonacci numbers in a Maple list L, step through a loop counting up the sequence, and in each loop iteration add to the list with L := [op(L), new_element].

f1 := proc(N)
    local i, L;
    L := [1, 1];
    for i from 2 to N do
        L := [op(L), L[-2] + L[-1]];
    end do;
    L
end proc:

The problem with this approach is that Maple lists are not intended to be grown in this way, and each redefinition of the list L costs O(N) time, making the whole procedure O(N).

> time(f1(N));
31.440

From now on we'll use only implementations that take O(N) time, and see how they compare.

The next implementation uses a Maple table stored as a local variable of a module. We compute the sequence and assign it to table entries, and at the end generate a Maple list from these table entries. The use of tables is advantageous because tables can be dynamically grown more efficiently than lists can, but this approach has the disdvantage that the table must be completely re-traversed in the last step to generate the list.

f2 := module()
    local ModuleApply, T;

    T := table([0=1, 1=1]);

    ModuleApply := proc(N)
        local i;
        for i from 2 to N do
            T[i] := T[i-2] + T[i-1];
        end do:
        [seq(T[i], i=0..N)]
    end proc:
end module:

You'll see that since this code avoids list-appending, it's quite a bit faster:

> time(f2(N));
1.110

The next idea is to do essentially what we did above, but using remember tables rather than built-in Maple tables. We'll pre-populate the remember table for the utility procedure p with entries for the base cases, F0 and F1, and then build our list through recursive calls.

f3 := module()
    local p, ModuleApply;

    p := proc(n) option remember;
        procname(n-2) + procname(n-1)
    end proc:

    p(0) := 1;
    p(1) := 1;
 
    ModuleApply := proc(N)
        [1, 1, seq(p(i), i=2..N)]
    end proc:
end module:

As we might expect, the computed result is pretty close (but a bit higher than) our last number.

> time(f3(N));
1.229

One disadvantage of the last two approaches is that they both rely on creating an intermediate data structure of size O(N) (the table or remember table), which is promptly discarded, since what we're interested in is a list. Instead we can use a 'cache', which is a new feature of Maple 9.5, and is essentially a remember table of bounded size.

f4 := module()
    local p, ModuleApply;

    p := proc(n) option cache;
        procname(n-2) + procname(n-1)
    end proc:

    p(0) := 1;
    p(1) := 1;

    ModuleApply := proc(N)
        [1, 1, seq(p(), i=2..N)]
    end proc:
end module:

It works out to be slightly slower than f3; this may be because of the time taken in garbage-collecting the information from the cache.

> time(f4(N));
1.459

Next, I've tried to implement something like a cache, without actually using one. Here we have a utility procedure with two counters, which represent the last and second-last Fibonacci number respectively. The procedure updates the values and returns the appropriate one. We'll build the list with recursive calls to this procedure. Like the cache, the intermediate data structures used here are of bounded size, since there are just two of them.

f5 := module()
    local p, n1, n2, ModuleApply;

    n1 := 1;
    n2 := 1;

    p := proc(n)
        (n1, n2) := (n2, n1+n2);
        n2;
    end proc:

    ModuleApply := proc(N)
        [1, 1, seq(p(), i=2..N)]
    end proc:
end module:

It's about the same speed as using the cache, which is impressive since it involves a lot more procedure calls into the procedure body than the cache example, which mostly consisted of calls into the cache.

> time(f5(N));
1.439

Finally, mostly for completeness, we try to see we get by pre-allocating an Array and adding to that. Arrays can be modified inplace in O(1) time, but cannot be dynamically grown. Here we make an array of the appropriate size, fill it, then convert it to a list.

f6 := proc(N)
    local A, i;
    A := Array(0..N, {0=1,1=1});
    for i from 2 to N do
        A[i] := A[i-2] + A[i-1];
    end do:
    convert(A, 'list')
end proc:

It's about as fast as the previous two examples, and has the advantage of simplicity of design. However, building the array and converting that to a list is another instance of using large intermediate data structures, which we avoided with the previous two examples.

> time(f6(N));
1.469

So it looks like f2 is the winner in this very specialized contest. This probably has more to do with the quirks of this example than any difference in efficiency of the code concerned. And if minimizing intermediate memory use is the goal, then f4 or f5 would be preferable.

First 19 20 21 Page 21 of 21