Advanced Search
All Posts : Product Tips & Techniques

The Eratosthenes Sieve

Alec Mihailovs 4475 Maple
September 08 2005
0
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, 
s := proc(n::integer[4], A::Array(datatype = integer[4]),
B::Array(datatype = integer[4]))::integer[4];
local k, j;
    for k from 2 to isqrt(n) do
        if A[k] = 0 then
            for j from k to iquo(n, k) do A[k*j] := 1 end do
        end if
    end do;
    j := 1;
    for k from 2 to n do
        if A[k] = 0 then B[j] := k; j := j + 1 end if
    end do
end proc:

cs := Compiler:-Compile(s):

ES := proc(n)
local A, B;
    A := Array(1 .. n, datatype = integer[4]);
    B := Array(1 .. numtheory:-pi(n), datatype = integer[4]);
    cs(n, A, B);
    B
end proc:
ES is much faster than f. For example, 
time(f(1000000));
                                13.438
time(ES(1000000));
                                0.155
The prime numbers produced by both procedures, are the same, 
is(convert(ES(1000000),set)=f(1000000));

                                 true
ES is also faster than the probabilistic prime procedure suggested by Stephen Forrest in his blog, 
sp5 := N->{seq(ithprime(i), i=1..numtheory:-pi(N))}:

time(sp5(1000000));

                                3.609

is(convert(ES(1000000),set)=sp5(1000000));

                                 true
time(ES(10000000));
                                1.671
time(sp5(10000000));
                               560.231
time(f(10000000));
                               505.578
As one may notice, f is also faster than s5 for large n. Another example, 
time(ES(2000000));
                                0.343
time(sp5(2000000));
                                44.312
time(f(2000000));
                                37.156
Even s is faster than s5 and f for large n, 
A := Array(1 .. 2000000, datatype = integer[4]):
B := Array(1 .. numtheory:-pi(2000000), datatype = integer[4]):
time(s(2000000,A,B));
                                14.623
By the way, ES works quite fast even for n=100,000,000. Looking at numbers above, I didn't try f or sp5 for that :) 
time(ES(100000000));
                                18.578
prime

Please Wait...