Erathostenes sieve (find all prime numbers lees than a given limit)

[Francesco Potorti`]

This routine has a functionality similar to that of find_primes, but
instead it finds all prime numbers not greater than n (find_primes
finds the first n prime numbers, excluding 1).

It uses the classical (a very suitable word !) Erathostenes sieve
algorithm, so often used as a quick and dirty compiler benchmark.

Its advantages are that it correctly includes 1 among the prime
numbers (find_primes forgets about that) and that it is a zillion
times faster (1 zillion is equal to about 3000 on an alpha for the
first 100000 primes).

Unfortuately, it also consumes much more memory, increasing linearly
with its argument.

-----------------------------------sieve.m-------------------------------
# Put in the public domain by:
# Francesco Potorti` <address@hidden> 1995

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

  if (!is_scalar(n) || n < 1 || n != round(n))
    error ("n should be a positive integer\n");
  endif
  if (n < 3)
    A = [1:n];
    return;
  endif
    
  N = [1, 2, 3:2:n];
  lN = length(N);
  limit = sqrt (n);
  prime = ones (size (N));
  i = 3;
  while (N(i) <= limit)
    prime(i+N(i):N(i):lN) = zeros (1, floor ((lN-i)/N(i)));
    while (N(i) <= limit && !prime(N(++i)))
    endwhile
  endwhile
  
  A = N (find (prime));
endfunction