Erathostenes sieve (even quicker)

[Francesco Potorti`]

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