\ factored sieve from
\ http://www.forthfreak.net/wiki/index.cgi?SieveOfEratosthenes

: sqrt ( u -- root )
    1 begin
	2dup / over + 2/ tuck
	= until
    nip ;

\ the 1000th prime is 7919 
\ 7919 2/ constant maxp      \ 2/ because we only count odd primes
8191 constant maxp

create flags maxp allot    \ one flag per odd number

: set-flags   flags maxp 1 fill ;

: count-flags ( -- n )
  0  flags maxp bounds do  I c@ +  loop ;

: clear-multiples ( prime start -- prime )
  flags maxp + swap do  0 I c!  dup +loop ;

: check-primes
  3 ( initial-prime )
  flags  maxp 2* sqrt 2/  bounds do
    I c@ if       \ prime?
      dup I + clear-multiples
    then 2 +      \ next odd number
  loop drop ;

: list-primes
  2 .      \ since we only count odd primes
  flags maxp bounds do  I c@ if  \ prime?
    I flags - 2* 3 + .
  then loop ;

: primes   set-flags check-primes ;
: count-primes ( -- n ) count-flags 1+ ;   \ +1 for 2

: benchmark ( -- )
    1000 0 ?do primes count-primes drop loop ;

\ calls
\  primes
\  count-primes .     \ 1000
\  list-primes        \ 2 3 5 7 11 .... 7919