-- demo/rosetta/Extensible_prime_generator.exw
sequence primes = {2,3,5,7}
atom sieved = 10

procedure add_block()
integer N = min((sieved-1)*sieved,400000)
sequence sieve = repeat(1,N)    -- sieve[i] is really i+sieved
    for i=2 to length(primes) do -- (evens filtered on output)
        atom p = primes[i], p2 = p*p
        if p2>sieved+N then exit end if
        if p2<sieved+1 then
            p2 += ceil((sieved+1-p2)/p)*p
        end if
        p2 -= sieved
        if and_bits(p2,1)=0 then p2 += p end if
--      if sieve[p2] then           -- dang!
            for k=p2 to N by p*2 do
                sieve[k] = 0
            end for
--      end if
    end for
    for i=1 to N by 2 do
        if sieve[i] then
            primes &= i+sieved
        end if
    end for
    sieved += N
end procedure

function is_prime(integer n)
    while sieved<n do
        add_block()
    end while
    return binary_search(n,primes)>0
end function

atom t0 = time()
while length(primes)<20 do add_block() end while
printf(1,"The first 20 primes are: ")   ?primes[1..20]
while sieved<150 do add_block() end while
sequence s = {}
for k=abs(binary_search(100,primes)) to length(primes) do
    integer p = primes[k]
    if p>150 then exit end if
    s &= p
end for
printf(1,"The primes between 100 and 150 are: ")    ?s
s = {}
for i=7700 to 8000 do
    if is_prime(i) then s&=i end if
end for
printf(1,"There are %d primes between 7700 and 8000.\n",length(s))
for i=1 to 8 do
    integer k = power(10,i)
    while length(primes)<k do
        add_block()
    end while
    printf(1,"The %,dth prime is : %d\n",{k,primes[k]})
end for
?time()-t0