[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
## 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.
--
Francesco Potorti` | address@hidden (Internet)
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

[Prev in Thread] |
**Current Thread** |
[Next in Thread] |

**Erathostenes sieve (even quicker)**,
*Francesco Potorti`* **<=**