
/* The *real* sieve of Erathosthenes. */

// Update an int vector in-place (no checks!).
using pointers;
infix 1000 := ;
def x!i := a = put_int (pointer x+i*SIZEOF_INT) a;

// The sieve. Returns the number of primes up to n.
sieve n = foldl sieve 0 (0..n-1)
// Uncomment this to get the vector of primes instead.
//	  $$ filter ((~=)0) x
with
  sieve count i = count+1 when
                    p = i+1; do (\j -> x!j := 0) (i+p:i+2*p..n-1);
                  end if x!i; // x!i is prime
                = count otherwise;
end when
  x = imatrix (1..n); x!0 := 0;
end;
