(       Title:  Eratosthenes Sieve
         File:  byteprimes.fs
   Derived by:  David N. Williams
      License:  Public Domain?
Starting date:  March 16, 2001
Version 0.7.1:  August 29, 2001
Version 0.7.3:  November 9, 2001
Version 0.7.4:  April 10, 2002
- Added test for finish > 2.

This implementation of the Eratosthenes sieve is essentially
that of J. Gilbreath, Byte Magazine, 9/81.  The version of the
original copied below is from the Forth, Inc. benchmark series.
In spite of their statement that it is the original, it is
probably modified slightly.
)

\ *** ORIGINAL VERSION

0 [IF]
\ =====================================================================
\
\ --------------------------------------------------------------------
\ Eratosthenes sieve benchmark program
\
\ This is the original BYTE benchmark.
\
\ --------------------------------------------------------------------

8190 CONSTANT SIZE
CREATE FLAGS   SIZE ALLOT

: DO-PRIME ( -- n )
   FLAGS SIZE -1 FILL
   0 SIZE 0 DO
      I FLAGS + C@ IF
         I 2* 3 + DUP I +  BEGIN
            DUP SIZE < WHILE
               DUP FLAGS + 0 SWAP C! OVER +
         REPEAT  2DROP  1+
   THEN  LOOP ;

: $SIEVE$ ( -- )
   BEGIN [$  DO-PRIME  SIZE $] UNTIL ." Eratosthenes sieve "
   .  ." Primes" ;
[THEN]

\ *** ERATOSTHENES SIEVE ALGORITHM
(
To find the primes from 3 up to and possibly including some odd
number, start with flags set for each of the odd numbers in the
range, and then successively clear those for multiples of each
prime in turn, starting with 3.  After each clearing sweep
starting at the flag for a prime, step from that flag to the
next one that remains set, which is guaranteed to correspond to
the next prime, and start the next sweep from there.

How many flags do we need?

If the input for the top of the range is an unsigned integer u,
the number of odd numbers <= u is u/2 if u is even or [u+1]/2 if
u is odd, i.e., [u+1] modulo 2.  Since we start with 3 instead
of 1, we actually need one less flag.

How far do we step from a prime multiple to get the next
multiple?

If we start from an odd multiple but include both even and odd
integers, the next two multiples [counting by the prime] will be
even and odd, respectively, two prime number steps to get to the
next odd multiple.  Since we only include flags for odd numbers,
we need half as many steps, or one prime number step in the list
of flags to get to the next one for an odd prime multiple.
)

\ *** PRINTING VERSIONS (ANS Forth)
(
There is no attempt here break the output into lines.  I use
these words mainly to remind me which numbers in a fairly small
range are primes, for various PIN's I've generated.
)

decimal

: .[primes]  ( u1 u2 -- )
(
Print the prime numbers in the inclusive range u1 <= u2, plus
the total number of primes printed.  All primes in the range 3
<= u2 are actually found, with only those in the required range
printed.
)
  ( start finish) 0 ( flags) 0 ( #flags) 0 ( limit)

  locals| limit #flags flags finish start |

  finish start u<
  ABORT" Range start must not be unsigned greater than end." 
  finish 3 u<
  ABORT" Range end must be larger than 2."

  finish 1+ 2/ 1- to #flags
  #flags allocate ABORT" Can't allocate prime flags."
      to flags
  flags #flags 1 fill
  flags #flags + to limit
  cr  \ start line of primes

  0                           ( #primes)
  limit flags DO
    i c@ IF  \ is prime
      i flags - 2* 3 +        ( _ prime)
      dup i +                 ( _ _ first.hole)
      BEGIN                   ( _ _ hole)
        dup limit u<
      WHILE
        ( hole) 0 over c!
        ( prime hole) over +  ( _ _ next.hole)
      REPEAT
      ( hole) drop
      ( prime) dup start u< IF ( prime) drop
          ELSE ( prime) . ( #primes) 1+ THEN
    THEN
  LOOP

  cr ( #primes) . ." primes "
  flags free ABORT" Error freeing prime flags." ;

: .primes  ( u -- )
(
Print the prime numbers from 3 up to possibly u, plus the total
number of primes printed.
)
  3 swap .[primes] ;