help-octave
[Top][All Lists]

## Erathostenes sieve (even quicker)

 From: Francesco Potorti` Subject: Erathostenes sieve (even quicker) Date: Thu, 03 Aug 1995 18:10 +0100 (MET)

```Another (the last one, I think) version of the Erathostene's sieve
algorithm fo rfinding all prime numbers less than a given integer.

This one is smaller, quicker and cleaner.  Moreover, Andreas
Weingessel made me note that indeed the ususal conventions say that 1
is not a prime number, so the output of sieve(100000) is equal to that
of list_primes(9592), but quicker by about 2700 times on an alpha.  It
also consumes much more memory.
--
researcher at          | 39369::pot (DECnet)
CNUCE-CNR, Pisa, Italy | +39-50-593203 (voice) 904052 (fax)

----------------------------------sieve.m---------------------------------
function A = sieve (n)
# Erathostenes sieve algorithm: find all primes not greater than n

# Put in the public domain by Francesco Potorti` <address@hidden> 1995

if (!is_scalar(n) || n < 1 || n != round(n))
error ("n should be a positive integer\n");
endif
if (n <= 3)
A = [2:n];
return;
endif

# Make an array of ones, one for the number 2 and one for each odd number
# greater than 2 and not greater than n.  The natural corresponding to the
# i-th one (i>2) is then 2*n-1.  Put a zero where the corresponding
# natural is not a prime.

lp = floor ((n+1)/2);
prime = ones (1, lp);
for i = 2:(sqrt(n)+1)/2
if prime(i)
Ni = 2*i-1;
range = (Ni*Ni+1)/2:Ni:lp;
prime(range) = zeros (size (range));
endif
endfor

A = 2 * find (prime) - 1;
A(1) = 2;
endfunction

```