Monday, November 23, 2015

Small Prime Sieve in FORTH

An explanation of the algorithm here can be found in a blog entry on the Sieve of Eratosthenes in my math and English blog, and in a previous entry in this blog on primes in C, along with my reasons for doing this, such as I have.

FORTH requires a lot of explanation, so I'm going to show you a C-like version first.

Comments in FORTH are between parentheses. Or you can start a comment with a backslash and it will end at the end of the line.

So, just read the code naturally and read the explanation in the comments. It's surprisingly readable (because I have written it to be readable). It's also full of hidden surprises for those used to algol-like in-fix languages including C and Ada.

Here's the C-like source:



( Archetypical implementation of the sieve of eratosthenes in FORTH -- gforth --
\ Copyright Joel Rees, 2015
\ By Joel Rees, Amagasaki, Japan, 2015
\ All rights reserved.
\ Permission granted by the author to use this code
\ for any purpose, 
\ on condition that substantial use 
\ shall retain this copyright and permission notice.
\
\ Timing results with MAXSIEVE set to 4194304 and output commented out:
\ me@fun:~/work/mathgames/sieve$ time gforth forthsieve.fs -e bye

\ real 0m1.061s
\ user 0m0.988s
\ sys 0m0.012s
\
\ Comparing to the C:
\
\ me@fun:~/work/mathgames/sieve$ time ./sieve_c
\ 
\ real    0m0.457s
\ user    0m0.436s
\ sys    0m0.020s
\
\ A bit more than double the run-times, 
\ so, really, in the same ballpark,
\ if not quite in the same league for speed.
)


256 constant MAXSIEVE
MAXSIEVE 1- 2 / constant FINALPASS

5 constant DISPWIDTH ( enough digits to display MAXSIEVE )

create sieve MAXSIEVE allot


: sieveMain ( -- )
0 sieve c!      ( 0 is not prime. )
0 sieve 1+ c!   ( 1 is not prime. )
sieve MAXSIEVE 2 ?do    ( set flags to true for 2 to FINALPASS. )
   -1 over i + c! loop  ( sieve pointer remains on stack. ) 
FINALPASS 2 ?do ( clear flags at all multiples. )
   dup i + c@ if      ( don't bother if not prime. )
      MAXSIEVE i dup + ?do    ( start at first multiple -- double. )
         0 over i + c!      ( clear at this multiple. )
      j +loop           ( sieve still on stack. )
   then
loop            ( sieve still on stack. )
MAXSIEVE 0 ?do
   i DISPWIDTH .r ." : is "
   dup i + c@ 0= if 
      ." not " 
   then
   ." prime." cr
loop drop ;

sieveMain




Okay, it looks like it ought to make sense, right?

But, if you are new to FORTH, and used to algol-like languages like C and Ada, it doesn't. So, let's look at it again, with comments that explain the FORTH itself.

But before we do that I should explain the stack.

FORTH's native interpreter operates in a post-fix grammar. In other words,
MAXSIEVE - 1
is written
MAXSIEVE 1 -
When the interpreter sees the constant MAXSIEVE, it pushes it on the stack. Then it sees 1 and pushes it on the stack. Then it sees the subtraction operator and subtracts the top item on the stack from the second:

inputstack
MAXSIEVEMAXSIEVE
1MAXSIEVE 1
-(result of MAXSIEVE - 1)

Or, borrowing from the code:
 MAXSIEVE 1- 2 / constant FINALPASS

inputstackwhat happens
MAXSIEVEMAXSIEVEConstant pushed on stack.
1-(result of decrement)Top item of stack gets decremented.
2(result1) 2 Numeric constant pushed on stack.
/(result of division)2nd on stack gets divided by top.
constant FINALPASS(empty)Constant FINALPASS generated in symbol table.

Now we can look at the code again:



256 constant MAXSIEVE
\ Defines MAXSIEVE as a true constant.
\ Incidentally, gforth is case insensitive, so we could call it "maxsieve", too.
\ But we won't.

MAXSIEVE 1- 2 / constant FINALPASS
\ See the explanation above.

5 constant DISPWIDTH ( enough columns to display MAXSIEVE )
\ We want this for integer output.

create sieve MAXSIEVE allot
\ This names and allocates a byte array named "sieve".
\ "sieve" in the input will now leave the array address on the stack:
\ ( -- adr )


: sieveMain ( --)
\ The colon operator defines a new word in the FORTH dictionary.
\ Or you can think of it as a new function or operator.
\ It is a symbol table entry with code that can be executed.
\ The ( -- ) is just a comment that shows tha there are no stack effects.
\ Some dialects of FORTH actually parse the stack effects, by the way.

\ "c!" stores the 2nd on stack at the address on top:
\ The stack effect is ( n adr -- )

0 sieve c!      ( 0 is not prime. )\ So, this stores a 0 in the first entry in the array,
\ for the integer zero.

0 sieve 1+ c!   ( 1 is not prime. ) \ And this stores a 0 in the second entry in the array,
\ for the integer one.


\ ?do 

\ Now we put the address on the stack again,
\ and set up a counted loop
\ with the limit and the initial count on stack
\ and start the loop, testing the limit against the count:
sieve MAXSIEVE 2 ?do    ( set flags to true for 2 to FINALPASS. )

\ "-1" is read as the numeric constant negative one.

\ "over" copies the 2nd on the stack and pushes it on top,
\ the stack effect is ( adr n -- adr n adr )

\ "i" pushes the current value of the loop counter on the stack.
\ The stack effect is ( -- counter )

\ "+" adds the top two numbers on stack, of course.
\ The stack effect is ( n1 n2 -- sum ).

\ "loop" marks the end of the loop.
\ It has no parameter stack effect: ( -- ),
\ but it removes the limit and current loop counter
\ from where they are stored.
\ It also branches back to the beginning of the loop,
\ if the limit has not been exceeded.

\ So, in the following code,
\ the loop counter is added to the array address,
\ and -1 is stored into the array as a byte value
\ at each byte in the array.

   -1 over i + c! loop  ( sieve pointer -- still on stack. )

\ This is the nested loop that clears the non-primes
\ by synthetic multiplication.
\ Note that "i" is always the innermost loop counter.
\ When we need to access the outer loop counter in the inner loop,
\ we do so with "j".

\ "?do" tests the limit against the initial value
\ and does not enter the loop if the limit is already exceeded.
\ "dup" duplicates the top of stack, effect is ( n -- n n ).
\ "c@" gets the byte value at the address on stack
\ and replaces the address with the value obtained.
\ Stack effect here is ( adr -- value )

\ "if" sets up a conditional evaluation.
\ Stack effect is ( flag -- ).
\ If the flag is zero, code up to "then" is skipped.
\ If the flag is non-zero, execution continues,
\ up to an optional "else".
\ The postfix-influenced interpretation of "then"
\ may be a little hard to get used to.
\ It is often aliased to "endif".

\ "+loop" is like loop, but it counts by the increment on top of stack,
\ stack effect is ( count -- )

FINALPASS 2 ?do ( clear flags at all multiples. )
\ Since the flag clearing pass starts at double each prime,
\ any prime over half of MAXSIEVE will never need to be looked at.

   dup i + c@ if      ( don't bother if not prime. )
      MAXSIEVE i dup + ?do    ( start at first multiple -- double. )
         0 over i + c!      ( clear at this multiple. )
      j +loop           ( sieve still on stack. )
\ Since j is the outer loop counter,
\ it is the current prime being scanned.
   then
loop            ( sieve still on stack. )

\ ".r" prints the 2nd integer on the stack
\ in a field width given by the top of stack.
\ "."" (sorry about the quotes) prints the following string,
\ up to the next quote.

\ "0=" inverts the flag on top of stack.
\ Non-zero become zero, zero becomes -1 .

\ "cr" writes a carriage return out.

\ "drop" drops a value from the stack,
\ stack effect is ( n -- ).

\ ";" terminates the colon definition.

MAXSIEVE 0 ?do
   i DISPWIDTH .r ." : is "
   dup i + c@ 0= if
      ." not "
   then
   ." prime." cr
loop drop ;

\ Now we can execute the newly defined function:

sieveMain
\ which generates the sieve and prints it.



As I said, this FORTH source looks like it was written by a C programmer. (It was. I am far more familiar with C than FORTH.)

You notice that there don't seem to be any local variables in FORTH. Modern FORTHs do have named local variables, but most local variables in FORTH are just parameters of expressions, anonymous on the stack.

Instead of local variables, FORTH traditionally cut functions up in smaller pieces. Some exercise of care allowed heavy re-use of those pieces, yielding more compact, highly readable code. (And then you could go to far and the code could become too dense to read.)

FORTH allows one to make new vocabularies, and that would allow us to hide the small pieces.

Finally, trying to use this on my bif-c revealed some of the bugs I haven't squashed yet. Trying to use it on the fig-FORTH for 6800 I transcribed some time back shows the wisdom of using the small pieces. The small pieces are far easier to debug. This works with the fig-FORTH, and, with the vocabulary stuff commented out, with my bif-c. (And now I have some debugging work to do on the bif-c.)


( Archetypical implementation of the sieve of eratosthenes in FORTH -- fig, bif-c --  )
( using more of the FORTH idiom. )
( Copyright Joel Rees, 2015 )
( By Joel Rees, Amagasaki, Japan, 2015 )
( All rights reserved. )
( Permission granted by the author to use this code )
( for any purpose,  )
( on condition that substantial use  )
( shall retain this copyright and permission notice. )



VOCABULARY sieve-local ( Make a local symbol table. )
sieve-local DEFINITIONS ( Switch to the local vocabulary. )


256 CONSTANT MAXSIEVE
MAXSIEVE 1 - 2 / CONSTANT FINALPASS

5 CONSTANT DISPWIDTH ( enough digits to display MAXSIEVE )


0 VARIABLE sieve ( Old FORTHs don't provide a default behavior for CREATE )
( gforth will leave the zero there. Old FORTHs need an initial value. )

   HERE sieve - DUP ( Old FORTHs don't provide a CELL width. )
   MAXSIEVE SWAP - ALLOT ( Allocate the rest of the byte array. )

   CONSTANT CELLWIDTH ( To show how this can be done. )

: sieveInit ( -- adr )
0 sieve C!      ( 0 is not prime. )
0 sieve 1+ C!   ( 1 is not prime. )
sieve MAXSIEVE 2 DO    ( set flags to true for 2 to FINALPASS. )
   -1 OVER I + C! LOOP  ( sieve pointer -- still on stack. )
;

: primePass ( adr prime -- adr )
MAXSIEVE OVER DUP + DO    ( start at first multiple -- double. )
   OVER I + 0 SWAP C!     ( clear at this multiple. )
   DUP +LOOP              ( next multiple )
DROP ;      ( sieve address still on stack. )

: findPrimes ( adr -- adr )
FINALPASS 2 DO   ( clear flags at all multiples. )
   DUP I + C@ IF ( don't bother if not prime. )
      I primePass
   ENDIF
LOOP ;           ( sieve still on stack. )

: printPrimes ( adr -- )
MAXSIEVE 0 DO
   I DISPWIDTH .R ." : is "
   DUP I + C@ 0= IF
      ." not "
   ENDIF
   ." prime." CR
LOOP DROP ;


FORTH DEFINITIONS

: sieveMain ( -- )
[ sieve-local ] sieveInit
findPrimes
printPrimes ;


sieveMain



If you are wondering about the speed of this version compared to the long version, it's pretty close. With the print commented out and the size set to 4194304, here's what we get in gforth:
me@fun:~/work/mathgames/sieve$ time gforth forthsieve.fs -e bye

real    0m1.096s
user    0m1.080s
sys    0m0.008s
which looks about right for 7 instructions vs. 6 in the inner loop.

(In bif-c, with the buggy parts commented out, too, it takes about 5 seconds. I guess, while I'm debugging it, I should add an option to parse an input file and go, so I can get accurate timings. But timings really aren't meaningful yet, there.)

You might also be interested in the FORTH language sieve entry in Rosetta Code.
Comparing the times for the current example using the dictionary buffer without allocating it:
me@fun:~/work/mathgames/sieve$ time gforth rs_sieve.fs -e bye

real    0m1.463s
user    0m1.444s
sys    0m0.012s
in gforth.

Friday, November 6, 2015

Small Prime Sieve in Ada (3 versions)

An explanation of the algorithm here can be found in a blog entry on the Sieve of Eratosthenes in my math and English blog, and in the previous blog entry, along with my reasons for doing this.

But that doesn't explain why I would do this three times.

Ada does not have a stepped loop. Altering the value of the index of a counted loop is considered dangerous, even if it is only to step the counter. Also, terminating a counted loop is considered difficult. So the counted loop is only present in forms that specify a range or a set of objects to count through. The programmer is to leave the details of setting up the count and counting through it to the compiler.

I understand the theory. I'll use strong language about this elsewhere, but I strongly disagree.

That said, there is one way to conform with the design of Ada. That is to set up a count (an artificial count in this case) and multiply by the step:


-- Archetypical implementation of the sieve of eratosthenes in Ada, with index multiply
-- By Joel Rees, Amagasaki, Japan, 2015
-- All rights reserved.
-- Permission granted by the author to use this code
-- for any purpose, 
-- on condition that substantial use 
-- shall retain this copyright and permission notice.
--
-- This is much faster than the modular division version, 
-- but reveals the effects of the tight semantics of Ada.
-- Still not anywhere close to the speed of the C version.
--
-- (With MAXSIEVE set to 131072:
-- me@fun:~/work/mathgames/sieve$ time ./adasieve >adaout.text
--
-- real    0m9.386s
-- user    0m1.112s
-- sys    0m7.948s
-- me@fun:~/work/mathgames/sieve$ time ./sieve_c > cout.text
--
-- real    0m0.148s
-- user    0m0.112s
-- sys    0m0.032s
--
-- With MAXSIEVE set to 1048576 and the output commented out:
-- me@fun:~/work/mathgames/sieve$ time ./adasieve 
-- 
-- real    0m0.125s
-- user    0m0.116s
-- sys    0m0.004s
-- me@fun:~/work/mathgames/sieve$ time ./sieve_c 
-- 
-- real    0m0.099s
-- user    0m0.096s
-- sys    0m0.004s
--
-- .. the latter of which is about the expected result.)
--
-- I read somewhere that it's the recommended approach.
--

with Ada.Text_IO; -- use Ada.Text_IO;    -- #include <stdio.h>
-- use would make it more like C, but I'm having trouble finding Integer_IO.
-- #include <stdlib.h> // We aren't using EXIT_SUCCESS

procedure AdaSieve is
    MAXSIEVE: constant := 256;
    MAXCANDIDATE: constant := MAXSIEVE - 1;
    FINALPASS: constant := ( ( MAXSIEVE + 1 ) / 2 );

    type Integer_Type is range 0 .. MAXCANDIDATE;
    subtype Candidate_Type is Integer_Type range
      2 .. Integer_Type'Last;

    package TextIO renames Ada.Text_IO;
    package IntIO is new Ada.Text_IO.Integer_IO( Integer_Type );
    sieve: array ( 2 .. MAXCANDIDATE ) of Boolean;
    -- index, prime: Integer_Type; -- loop variables
    -- This is the recommended approach!
    limit: Integer_Type;

begin
    -- sieve[ 0 ] := 0;
    -- sieve[ 1 ] := 0;
    for index in sieve'Range loop
        sieve( index ) := True;
    end loop;

    for prime in 2 .. FINALPASS loop
        if sieve( prime ) = True then
            limit := Integer_Type( Standard.Integer( Integer_Type'last ) / prime ); -- (gulp!)
            -- index check failed:
            -- limit := Integer_Type( ( Standard.Integer( Integer_Type'last ) + 1 ) / prime );
            for index in 2 .. Integer( limit ) loop 
-- TextIO.Put( "clearing index " ); -- debug
-- IntIO.Put( Item  => Integer_Type( index ) ); -- debug
-- TextIO.Put( " element " ); -- debug
-- IntIO.Put( Item  => Integer_Type( index * prime ) ); -- debug
-- TextIO.Put_Line( "." ); -- debug
                sieve( index * prime ) := False; 
            end loop;
        end if;
    end loop;

    
    IntIO.Put( Item  => 0 );
    TextIO.Put_Line( ": is not prime." );
    IntIO.Put( Item  => 1 );
    TextIO.Put_Line( ": is not prime." );
    for index in sieve'Range loop
        IntIO.Put( Item  => Integer_Type( index ),
                   Width => IntIO.Default_Width,
                   Base  => IntIO.Default_Base );
        TextIO.Put( ": is " );
        if not sieve( index ) then
            TextIO.Put( "not " );
        end if;
        TextIO.Put_Line( "prime." );
    end loop;
    -- return EXIT_SUCCESS;
end AdaSieve;



This approach is useful for walking through the bitbuckets in certain hash implementations, but I don't consider it a natural fit for this algorithm. There is no particular motivation for the division that tells how many times the loop must run, and the multiplication, while much less time consuming than division, is not as cheap as the add.

(You'll note the debugging prints I used to make sure the count was surviving the typecasting.)


Then there is a way to cheat and use the natural range, but use modular division to skip the steps where the loop should not execute. You get punished for using this one by rather slow running time:


-- Archetypical implementation of the sieve of eratosthenes in Ada, with index modular divide
-- By Joel Rees, Amagasaki, Japan, 2015
-- All rights reserved.
-- Permission granted by the author to use this code
-- for any purpose, 
-- on condition that substantial use 
-- shall retain this copyright and permission notice.
--
-- The modular division in the inner loop makes this really, really slow.
--

with Ada.Text_IO; -- use Ada.Text_IO;    -- #include <stdio.h>
-- use would make it more like C, but I'm having trouble finding Integer_IO.
-- #include <stdlib.h> // We aren't using EXIT_SUCCESS

procedure AdaSieve is
    MAXSIEVE: constant := 256;
    MAXCANDIDATE: constant := MAXSIEVE - 1;
    FINALPASS: constant := ( ( MAXSIEVE + 1 ) / 2 );

    type Integer_Type is range 0 .. MAXCANDIDATE;
    subtype Candidate_Type is Integer_Type range
      2 .. Integer_Type'Last;

    package TextIO renames Ada.Text_IO;
    package IntIO is new Ada.Text_IO.Integer_IO( Integer_Type );
    sieve: array ( 2 .. MAXCANDIDATE ) of Boolean;
    -- index, prime: Integer_Type; -- loop variables

begin
    -- sieve[ 0 ] := 0;
    -- sieve[ 1 ] := 0;
    for index in sieve'Range loop
        sieve( index ) := True;
    end loop;

    for prime in 2 .. FINALPASS loop
        if sieve( prime ) = True then
            for index in prime + prime .. MAXCANDIDATE loop 
                if index mod prime = 0 then
                    sieve( index ) := False;
                  end if;
            end loop;
        end if;
    end loop;

    
    IntIO.Put( Item  => 0 );
    TextIO.Put_Line( ": is not prime." );
    IntIO.Put( Item  => 1 );
    TextIO.Put_Line( ": is not prime." );
    for index in sieve'Range loop
        IntIO.Put( Item  => Integer_Type( index ),
                   Width => IntIO.Default_Width,
                   Base  => IntIO.Default_Base );
        TextIO.Put( ": is " );
        if not sieve( index ) then
            TextIO.Put( "not " );
        end if;
        TextIO.Put_Line( "prime." );
    end loop;
    -- return EXIT_SUCCESS;
end AdaSieve;



This is actually a useful pattern. There are numerous algorithms in which one examines every member of a set or array in turn, performing some calculation to determine how to process the element under consideration.


Then there is a work-around, manufacturing the counted loop using a while. Yes, if you can handle counted loop termination, you also know how to manufacture a counted loop with a general pre-test loop:


-- Archetypical implementation of the sieve of eratosthenes in Ada, with while loop
-- By Joel Rees, Amagasaki, Japan, 2015
-- All rights reserved.
-- Permission granted by the author to use this code
-- for any purpose, 
-- on condition that substantial use 
-- shall retain this copyright and permission notice.
--
-- Timing results with MAXSIEVE set to 4194304 and output commented out:
-- me@fun:~/work/mathgames/sieve$ time ./adasieve
--
-- real    0m0.464s
-- user    0m0.448s
-- sys    0m0.012s
-- me@fun:~/work/mathgames/sieve$ time ./sieve_c
--
-- real    0m0.457s
-- user    0m0.436s
-- sys    0m0.020s
--
-- This is what we should expect for a running time.
--

with Ada.Text_IO; -- use Ada.Text_IO;    -- #include <stdio.h>
-- use would make it more like C, but I'm having trouble finding Integer_IO.
-- #include <stdlib.h> // We aren't using EXIT_SUCCESS

procedure AdaSieve is
    MAXSIEVE: constant := 256;
    MAXCANDIDATE: constant := MAXSIEVE - 1;
    FINALPASS: constant := ( ( MAXSIEVE + 1 ) / 2 );

    type Integer_Type is range 0 .. MAXCANDIDATE;
    subtype Candidate_Type is Integer_Type range
      2 .. Integer_Type'Last;

    package TextIO renames Ada.Text_IO;
    package IntIO is new Ada.Text_IO.Integer_IO( Integer_Type );
    sieve: array ( 2 .. MAXCANDIDATE ) of Boolean;
    windex: Integer; -- Not Candidate_Type, to use it as array index?

begin
    -- sieve[ 0 ] := 0;
    -- sieve[ 1 ] := 0;
    for index in sieve'Range loop
        sieve( index ) := True;
    end loop;

    for prime in 2 .. FINALPASS loop
        if sieve( prime ) = True then
            -- Manufacturing the for loop with while
            windex := prime + prime;
            while windex < MAXSIEVE loop 
                sieve( windex ) := False;
                windex := windex + prime;
            end loop;
        end if;
    end loop;

    
    IntIO.Put( Item  => 0 );
    TextIO.Put_Line( ": is not prime." );
    IntIO.Put( Item  => 1 );
    TextIO.Put_Line( ": is not prime." );
    for index in sieve'Range loop
        IntIO.Put( Item  => Integer_Type( index ),
                   Width => IntIO.Default_Width,
                   Base  => IntIO.Default_Base );
        TextIO.Put( ": is " );
        if not sieve( index ) then
            TextIO.Put( "not " );
        end if;
        TextIO.Put_Line( "prime." );
    end loop;
    -- return EXIT_SUCCESS;
end AdaSieve;



Even though the loop is a while loop instead of a for loop, this most closely matches my original C source.

The FORTH implementation can be found here.

[You might also be interested in the Ada language sieve entry in Rosetta Code.]

Tuesday, November 3, 2015

Small Prime Sieve in C

You can find primes without using multiplication or division. All you need is a large place to write a large table of numbers.
I showed how to do the brute force method and the easy sieve of Eratosthenes by hand in my Math and English blog, here.

The easy version of the sieve of Eratosthenes is pretty simple. It's practically de rigueur. If you are a computer professional and you haven't yet written one out from scratch, you should look at yourself in the mirror and then sit down and do it now.

If you are not a computer professional, you can copy the source codes I'm going to present here into a source code editor and try them. Maybe it will get you interested in math and programming, or maybe you'll just say, "huh?" But it will be a better way to waste an hour or two than getting stoned.

My motivation in writing these is that I need to re-learn FORTH, and I think I should learn some Ada, and this is one of those problems that, when done the interesting way, exposes a lot of the basic features of a language.

This post is not going to expose much, however. It's just going to get us started.

The algorithm, to a fair degree of detail, will be as follows:

Set up a sieve array. Note that, in the representation below, the numbers are not in the array itself. They are indexes, or addresses into the array. The contents are the lower half.

In C:

    char sieve[ 24 ];

gives us an array called sieve with indexes from 0 to 23:

01234567 89101112131415 1617181920212223
________ ________ ________

By the underscores ("_") here, I mean that we don't know what is in the array. So the first step is to set the sieve array to a state that we can use.

We'll flag non-primes with "0".

Zero and one are by definition, non-prime. (They are integers, and they are not evenly divisible by any other integer, but they are not greater than 1.)

All the rest, we will say, may be prime. We don't know yet. Or, at least, the computer doesn't know yet. (Computers really aren't smart until we program some smarts into them.) So we will set their sieve entries to "1" so we can program some smarts into your computer.

(And it will forget, as soon as the program ends, unless we tell it to write it down somewhere on a disk drive or in flash memory or such. That's the way computers are.)

In C:

    sieve[ 0 ] = 0;
    sieve[ 1 ] = 0;
    for ( index = 2; index < 24; index = index + 1 )
    {   sieve[ index ] = 1;
    }

01234567 89101112131415 1617181920212223
00111111 11111111 11111111

The second step is to decide we are going to start at two.

    prime = 2;

The third step is to start at the first multiple of two and count up through the array, clearing the array at every multiple:

    4, 6, 8, 10, 12, 14, 16, 18, 20, 22

In C:

    for ( index = prime + prime; index < 24; index = index + prime )
    {   sieve[ index ] = 0;
    }

The result of this step when prime = 2 looks like this:

01234567 89101112131415 1617181920212223
00110101 01010101 01010101

We start at 2 times 2, then repeatedly add 2 to the index and clear the sieve array at the resulting multiple of 2.

It will be the same at 3:

    prime = prime + 1;
    for ( index = prime + prime; index < 24; index = index + prime )
    {   sieve[ index ] = 0;
    }

After the pass for 3, the array will look like this:

01234567 89101112131415 1617181920212223
00110101 00010100 01010001

By using prime to count through the potential primes, and index to clear the multiples, we can use the same loop for 2, 3, and all the primes we try.

Really, really simplifying, we will find that every multiple of a non-prime is also a multiple of a prime, so we could just count the variable we call "prime" through the whole list.

    for( prime = 2; prime < 24; prime = prime + 1)
    {
        for ( index = prime + prime; index < 24; index = index + prime )
        {   sieve[ index ] = 0;
        }
    }

When we are testing any prime that is greater than or equal to 12, the first multiple will be greater than or equal to 12 times 2, or 24, which is beyond the end of the array, outside our list.

Fortunately, the inner loop that actually clears the array will check whether index less than 24 before it enters the body of the loop the first time. So the loop won't even start once unless prime is less than 12.

We don't have to in this version, and it doesn't really save much time (relatively speaking) but we can go ahead and stop the outer loop at that point:

    for( prime = 2; prime < 12; prime = prime + 1)
    {
        for ( index = prime + prime; index < 24; index = index + prime )
        {   sieve[ index ] = 0;
        }
    }

If we knew that the flag in the sieve were always valid for the value of prime, then we could skip the inner loop when the flag is "0". This could potentially save a lot more time than just stopping the outer loop early.

Let's walk through the loops a little more carefully.

When prime is 2, the inner loop will clear all multiples of 2, In the process, it will also clear all multiples of multiples of 2:
  • 4 = 2 × 2
  • 6 = 2 × 3
  • 8 = 2 × 4 = (2 × 2) × 2
  • 10 = 2 × 5
  • 12 = 2 × 6 = (2 × 3) × 2
  • 14 = 2 × 7
  • 16 = 2 × 8 = (2 × 4) × 2
  • 18 = 2 × 9 = (3 × 3) × 2
  • 20 = 2 × 10 = (2 × 5) × 2
  • 22 = 2 × 11
So, when the pass for 2 is finished, effectively, so are the passes for 4, 6, 8, and 10. None of those passes have to actually be run.

Likewise, when the pass for 3 has finished, the passes for 6 and 9 don't have to be run. And when the pass for 5 has finished, the pass for 10 doesn't have to be run. How can we take advantage of this information?

After initially setting up the sieve, the flags for 2 and 3 are valid. So, when prime is 2, the flag is valid. The inner loop won't touch 3 when prime is 2, so when prime is 3, the flag is valid.

The inner loop clears the flag at 4 when prime is 2, so, by the time prime is 4, the flag is valid.

Looking carefully, we can see that, for every value of prime, the flags between prime and prime times 2 are always valid. The computer doesn't know, but we do, so we can tell the computer to use that information, by preventing the inner loop when the flag says it isn't really prime:

    for( prime = 2; prime < 12; prime = prime + 1)
    {
        if ( sieve[ prime ] == 1 )
        {   for ( index = prime + prime; index < 24; index = index + prime )
            {   sieve[ index ] = 0;
            }
        }
    }

So we now have a nice framework to work within.

Here is the source for a C version able to find all primes less than 256, with all the variables declared, some convenient manifest constants to make the size of the list easier to change, and other stuff I usually do when I write in C:



/* Archetypical implementation of the sieve of eratosthenes in C
** By Joel Rees, Amagasaki, Japan, 2015
** All rights reserved.
** Permission granted by the author to use this code
** for any purpose, 
** on condition that substantial use 
** shall retain this copyright and permission notice.
*/

#include <stdio.h>
#include <stdlib.h>


#define MAXSIEVE 256

#define FINALPASS ( ( MAXSIEVE + 1 ) / 2 )


int main( int argc, char * argv[] )
{   char sieve[ MAXSIEVE ];
    int index, prime;

    sieve[ 0 ] = 0;
    sieve[ 1 ] = 0;
    for ( index = 2; index < MAXSIEVE; index = index + 1 )
    {   sieve[ index ] = 1;
    }

    for ( prime = 2; prime < FINALPASS; prime = prime + 1)
    {
        if ( sieve[ prime ] == 1 )
        {   for ( index = prime + prime; index < MAXSIEVE; index = index + prime )
            {   sieve[ index ] = 0;
            }
        }
    }

    for ( index = 0; index < MAXSIEVE; index = index + 1 )
    {   printf( "%d: %s prime.\n", index, ( sieve[ index ] == 0 ) ? "is not" : "is" );
    }
    return EXIT_SUCCESS;
}



This is a long post, so my Ada and FORTH implementations will be posted separately:
  • Three versions of the sieve in Ada
    Showing three minor variations of the above algorithm,
  • The sieve in FORTH, but not natural FORTH. Includes code that will compile in 16 kilobytes of RAM on old 8-bit microcomputers, but, of course, won't give you large sieves.

[You might also be interested in the C language sieve entry in Rosetta Code.]