( 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] ;