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