(       Title:  Eratosthenes Sieve
         File:  primes.fs
   Derived by:  David N. Williams
      License:  Public Domain [if the original code is]
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.
)

\ *** ORIGINAL VERSION
(
This nonprinting version is from the Gforth benchmark file
sieve.fs.  It finds 1899 primes in the range 3 through 16381.
)

0 [IF]
DECIMAL
: SECS TIME&DATE  2DROP DROP  60 * + 60 * + ;
CREATE FLAGS 8190 ALLOT
FLAGS 8190 + CONSTANT EFLAG

: PRIMES  ( -- n )  FLAGS 8190 1 FILL  0 3  EFLAG FLAGS
  DO   I C@
       IF  DUP I + DUP EFLAG <
           IF    EFLAG SWAP
                 DO  0 I C! DUP  +LOOP
           ELSE  DROP  THEN  SWAP 1+ SWAP
           THEN  2 +
       LOOP  DROP ;
[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) 0 ( #primes) 3 ( first.odd)

  locals| odd #primes flags limit #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 + to limit
  flags #flags 1 fill

  cr  \ start line of primes
  limit flags DO
    i c@ IF  \ odd is prime
      i odd + ( hole) limit u< IF
        limit i odd +
        ( limit hole) DO 0 i c! odd +LOOP
      THEN
      odd start u< 0= IF
        odd ( prime) .   #primes 1+ to #primes
      THEN
    THEN
    odd 2 + to odd
  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] ;