Eratosthenes lived about 2200 years ago. One of his achievement in Mathematics was his method of generating prime numbers, called the Sieve of Eratosthenes. This is a very elegant way of quickly identifying which numbers are prime. It is also very easily and simply made into a computer program.

I have been experimenting with this method, and am absolutely stunned at how fast it can be. My first test was calculating the primes between 1 and 1,000,000. This executed so quickly, it took me 30 minutes of investigation to realize that it really only took 0.16 seconds to execute! Since then I have reduced the execution down to the point where the windows overhead within the program overshadows the actual calculation time, which became un-measurable on such a small number of primes. Recent results on my PIII with 394 megs of ram include

1-10,000,000 in under 1 second,

1-100,000,000 in about 10 seconds

1-1,000,000,000 in 122 seconds

and 1-5,000,000,000 is 497seconds.

The program is written using Delphi 7, and the coding language is Pascal. My first attempt was a classic implementation, with no shortcuts. It works on the assumption that 1 is not prime. An array of Boolean fields is created, one for each number in the range. A first loop loops from 1 to the square root of the limit of the range, and when it finds a prime, sets all multiples of that prime to 'not prime', by incrementing the position by the prime value until the upper limit is passed. This continues until all non primes are identified. The Pascal code is reproduced below. Feel free to criticise!

Procedure sieveold; begin {sets length of array to variable from form, and factor limit to square root of limit} SetLength(bitarray, arraylength); sqrtlength := round(sqrt(arraylength)); itarray[1] := 1; for i := 1 to sqrtlength do if bitarray[i] = 0 then begin j := i*2; while j <= arraylength do begin bitarray[j] := 1; j:= j + i; end; end; end;

The next version did a little bit of optimisation, including setting all even numbers to non prime in a separate pass, and the setting the multiples of the odd primes missing out all the even multiples. This reduced the operating time by 50% or so. Example code follows.

procedure sieve; var i:integer; begin SetLength(bitarray, arraylength); sqrtlength := round(sqrt(arraylength)); bitarray[1] := 1; {removes even numbers > 2 from array} j:=4; while j <= arraylength do begin bitarray[j] := 1; j:= j + 2; end; {for each prime (byte not true) set all subsequent odd multiples of it to true} for i := 1 to sqrtlength do if bitarray[i] = 0 then begin j := i*i; while j <= arraylength do begin bitarray[j] := 1; j:= j + i + i; end; end; end; end;

The main limitation with this version appeared to be the maximum size, which is restricted to the amount of available memory If your limits exceed the maximum amount (over 250,000,000 on my machine), it starts to swap memory in and out to disk. As you are hitting the WHOLE TABLE at least 15,000 times when the table gets this large, it tries to do a lot of swapping!

As I wanted to get up to 1,000,000,000, I decided to change the table to an array of bits rather than byte booleans. Pascal has the usual bitwise processing available, so I decided I would sacrifice some speed to get the bigger range (2,000,000,000 would be possible). To my surprise, swapping one instruction that set a Boolean field to true with 4 logical statements to divide the position by 8, giving the byte number and position with the byte, and to create a mask for the correct position, and ANDing or ORing it with the appropriate byte, was actually quicker, halving the execution time, as well as extending the range 8 fold! Coding follows.

procedure setbit (pos:integer); begin bytepos := pos shr 3; bitpos := pos and mask7; bitmask := mask1 shl bitpos; bitarray[bytepos] := bitarray[bytepos] or bitmask; end; function checkbit (pos:integer):bool; begin bytepos := pos shr 3; bitpos := pos and mask7; bitmask := mask1 shl bitpos; checkbit := BOOL(bitarray[bytepos] and bitmask); end; procedure sieve; {calculation routine separated out, only this part is timed} var i:integer; begin {sets length of array to variable from form, and factor limit to square root of limit}

SetLength(bitarray, (arraylength div 8)+1); sqrtlength := trunc(sqrt(arraylength)); mask7 := 7; mask1 := 1; setbit(1); {removes even numbers > 2 from array} j:=4; while j <= arraylength do begin setbit(j); j:= j + 2; end; {for each prime (byte not true) set all subsequent odd multiples of it to true} i:=3; while i < sqrtlength do begin if not checkbit(i) then begin j := i*i; while j <= arraylength do begin setbit(j); j:= j + i + i; end; end; inc(i,2); end; end;

And next, a version that goes as high as I can while remaining in memory, rather than using disk storage. This version acknowledges that a. 1 is a not a prime number, b. that 2 is a special case, and does not need to be included in the general case, and that all even numbers after 2 are not prime (without having to generate that information). This allows me to put 16 numbers into each byte in effect, so my spare 250,000,000 bytes of memory can now be used to calculate primes up to 4,000,000,000. The bitwise logic becomes a bit more complicated, dividing the position by 16 not 8 to give the correct byte, but dividing the remainder (position mod 16) by 2 to give the bit position. An extra SHR is used to get the bitpos. As we have broken to 2 billion barrier (well 2**31 anyway) it became necessary to change all integer variables from integer 32 to integer 64. Not a big deal, but the SQRT function does not work on int64, so was changed to a POWER function, using 0.5 as the exponent. It's probably faster anyway!

So, this version is just as fast for the lower numbers, with the added advantage of calculating the 189,961,812 prime numbers less than 4,000,000,000 is just a touch under 8 minutes (478.459 seconds).

The code looks a bit like this

procedure setbit (pos:int64); begin bytepos := pos shr 4; bitpos := pos and mask15; bitpos := bitpos shr 1; bitmask := mask1 shl bitpos; bitarray[bytepos] := bitarray[bytepos] or bitmask; end; function checkbit (pos:int64):bool; begin bytepos := pos shr 4; bitpos := pos and mask15; bitpos := bitpos shr 1; bitmask := mask1 shl bitpos; checkbit := BOOL(bitarray[bytepos] and bitmask); end; procedure TCalculatePrimes.GoClick(Sender: TObject); var i,j,sqrtlength : int64; code: integer; stri, strj, strl:string; starttime,endtime:tdatetime; elapsed: real; procedure sieve; {calculation routine separated out, only this part is timed} var i:int64; begin {sets length of array to variable from form, and factor limit to square root of limit} SetLength(bitarray, (arraylength div 16)+1); sqrtlength := trunc(power(arraylength, 0.5)); mask15 := 15; mask1 := 1; {for each prime (byte not true) set all subsequent odd multiples of it to true} i:=3; while i < sqrtlength do begin if not checkbit(i) then begin j := i*i; while j <= arraylength do begin setbit(j); j:= j + i + i; end; end; inc(i,2); end; end;

A colleague at work suggested that embedding the main "setbit" routine into the main code woould reduce the run time by 30-60%, which I seriously doubted, but I had a go anyway. The improvent in speed was about 10%, worthwhile, but not that dramatic. I tried putting all the bitwise logic statements into a single statement, but that slowed it down again!

The last nice touch was to calculate the amount of memory available, and to determine the highest number which could be calculated without the application going outside the memory, causing lost of disk activity. Once again I was surprised how much memory was available in my machine with windows running.

I managed to get a run with 5,000,000,000 as the upper limit, and calculated all the 234,954,223 primes between 1 and 5,000,000,000 in 497.019 seconds! A great programe for pushing a PC to its limit. I can raise the CPU temperature from 29 to 36 degrees while it is running!

Most of the above was written in 2004, and the application was then written in Borland Delphi. Borland Delphi is no more, but exists in another form as Embarcadero Delphi. I tried my very hardest, but could not manage to download and try a version of this software, despite being willing to possibly buy a copy if it turned out to be OK. So, my clever daughter found Lazarus, an open source Delphi clone which uses Free Pascal as its base.

I have converted the project from Borland Delphi to Lazarus (or rather lazarus did it for me with a couple of clicks and two small corrections to my code). I took the opportunity to improve the code a bit, and you can download the complete Lazarus project here and a copy of the Lazarus development tool here.

The latest version, run on a more modern machine with more ram is a bit faster, calculating all 455,052,511 primes between 1 and 9,999,999,999 (10 billion) in just 231 seconds (under 4 minutes). The screen print above is of the new version.

Fancy trying it? Windows users download Primes below, unzip it and execute it. It shows the time taken, number of primes in the range, and the highest prime, and gives you two options for asking about primes. Put a number in either field, and tab out, and it will give you the results.

Hopefully the next step will be an android version, which is a whole new area of learning for me, which Lazarus should help a lot with.

Dave Glover, 27/01/2014

Other stuff I may want to download!

CSV