How to sieve prime numbers?

Spoken mathematically, (for all numbers n):((for all numbers i>n): n mod i <> 0). And, the definition of a prime number is: (n is prime number) <=> ((for all numbers 1<i<n):n mod i <> 0).

Now imagine an array of all numbers from 2 to any upper limit, say, 10 million. Each number has a flag whether the number is checked or not. The array is initialized with all numbers unchecked.

You now begin with two and proceed upwards. Whenever you encounter an unchecked number, it is a prime number. The trick is, you now check out all multiples of this number, that is, all numbers that are dividable by this number and which are therefore not prime.

Now whenever you get to a number n, you have checked out all multiples of 1<i<n, and if n is not checked, it is therefore prime. You'll understand the algorithm best if you take a piece of paper and do it by hand maybe for 1..100.

Of course, there are two problems:

1. You need an array for the flags of all numbers in the interval from which you want to extract the primes. If you want to calculate all primes from 1 to 1.000.000, you need an array with 1.000.000 elements.
2. You must begin the algorithm with 2. The interval must not necessarily start with 1, but nevertheless you must begin the ruling-out with all potential dividents.

The below program does implement this sieve, plus a few goodies and optimizations.

• It doesn't check even numbers. The only even prime number is 2. And we don't need to rule out multiples of even numbers, because multiples of even numbers are even and therefore cannot be prime.
• Only one bit is allocated as a flag for each number. Together with the above optimization, each bytes holds information for 16 numbers. To calculate 1000 additional primes, 63 additional bytes of memory are needed.
• The program was originally developed on DOS, and had to suffer from the 640kB barrier. To minimize suffering, it checks upon startup which array size it is allowed to allocate, and adjusts itself appropriately. It does this by trying to allocate a much too big chunk of memory (for 10.000.000 numbers, needing about 640k) and then making the request smaller step by step until the allocation is granted.
• It first calculates all the primes, and when it is finished, it will prompt you for a display interval, and then it displays all primes within this interval.

I have to admit, because of not handling some special cases, it will identify 2 as not prime and 1 as prime. While it is clearly wrong about 2, we could discuss about 1 being prime or not ...

Compiling and Running

The program takes an upper limit as parameter, i.e. it will compute prime numbers up to this number given on the command line. This upper limit defaults, for historic reasons, to 14.000.000.

This should take substantially less than one minute on a workstation and on fast PCs (our HP/UX's here need 23 seconds). It will then prompt you for the beginning and end of an interval (you may enter even numbers), display all prime numbers in this area and then prompt again. The program terminates when you enter an end of interval that is less or equal the beginning of the interval.