Bounded Buffer Sieve of Erastophenes

• Download source - 2.7 KB

Introduction

When a stream of primes is required, typically, programmers use the Sieve of Erastophenes[1] algorithm to generate the numbers. A traditional sieve works by allocating a buffer the size of the largest prime being sought. For example, if one needs primes up to n, then the buffer is allocated from 0 to n-1. The sieve then works by starting with the number two and marking every multiple of two in the buffer as being composite. The output is then the number two. For the next iteration, it finds the next and first number not marked as composite, and starts setting all of its multiples as composite. It then reports this number, which is three. The algorithm proceeds this way until no more numbers can be found which are not marked composite up to n.

Background

The central problem with the traditional sieve is that the buffer needed is the size of the maximum prime required. For example, if 1 million is the largest prime required, then the buffer will be 1 million elements large and half of all the values in the table will be composite. To remedy this issue, a new sieve was designed which uses a fixed size buffer. The buffer size can be less than the largest prime desired. The algorithm works by collecting the primes that are discovered in a linked list. Once the buffer is exhausted, it resets the buffer, then uses the primes in the linked list to mark the composites in the buffer. When the buffer is exhausted, it repeats this process. The algorithm correctness comes from the fact that for any composite integer, there must exist a prime number that is less than or equal to the ceiling of the composite integer's square root. The prime numbers collected in the linked list suffice this requirement for the next set of composites in the buffer. This allows for discovery of prime numbers without necessarily allocating an n sized buffer. Since the number of primes is much less than the number of composites, the sieve can generate a larger number of primes than with the classic sieve.

The new algorithm can be divided into two sections. The first section is based on a traditional sieve to generate the first set of prime numbers which are collected into a linked list. The second section is the generation of primes using the discovered primes list. The key to understanding the algorithm is the following loops:

prime_tab[0] = true;
register integer lstart = (UNIT * (next-ONE));
register integer index  = 0;
for(list_int_t::iterator i = prime_list->begin();
                       i != prime_list->end();
                       ++i ) {
       	/* Mark which ones in unit * next
       	 * space are multiples of known primes.
       	 */
       	for(register integer j = ((lstart/(*i)) + ONE) * (*i);
                      (j - lstart) < UNIT; j += (*i) ) {
               	index = j - lstart;
               	prime_tab[index] = true;
       	}
}

The linked list of primes is iterated in the outer for loop. Then in the inner for loop, a given prime number is used to determine a starting index within the buffer (prime_tab) which is the starting integer to be marked as composite. Thereafter, all multiples of that prime within the table are marked as composite. The integer UNIT refers to the size of the buffer and next is the number of passes through the buffer. Since we are reusing the buffer, the index representing a slot in the buffer must logically map to the actual composite it represents. Therefore, if the buffer was four elements large (UNIT==4), then a pass means we have finished iterating it, marked all the numbers that are composite, and discovered all the primes. On the next pass, we then clear out our settings (reset prime_tab), and repeat. Except that the index 0 at that point would not represent integer 2, but would represent integer 6 as the first pass through the buffer considered integers 2, 3, 4, and 5.

In the code above, prime_tab[0] is always composite given that UNIT is designed to be divisible by two and is always initialized to true (composite). lstart is the logical start and this is computed as (UNIT * (next-ONE)) which is the logical value that the zero slot in the prime_tab buffer represents. In the inner for loop, integer j represents a logical slot and is initialized using integer division by dividing the logical start by the next prime in the list (*i), adding one to it (which has the effect of taking the ceiling of the result of dividing the logical start by the prime) and then multiplying this factor by the prime to get the logical slot. The actual slot is j - lstart which must be less than the size of the prime_tab buffer UNIT. Incrementing j by the prime (*i) at each iteration allows the algorithm to mark out all multiples of the prime for that logical pass. Note that for UNIT * next we compute lstart as (UNIT * (next-ONE)), with the algorithm starting with the number of passes (next) equal to 1. Considering our example of a buffer where UNIT==4, and pass (next) is 2, we compute: (lstart = (4 * (2-1)))==4. There would be three known primes in the prime_list (2,3,5). Using prime number 2 as an example, we have:(j = (4/2)+1 * 2)==6 and that 6 is the logical slot that 0 represents in the second pass. The physical slot is j - lstart which is 6-4==2. Thus actual slot 2 would be marked as composite. Then j would be incremented by the prime 2 to 8. The next actual slot would be 8-4==4. Since this is not less than UNIT==4, the algorithm would stop processing for that prime. It is possible that UNIT can be too small and the algorithm was tested with UNIT==1000 as the smallest value.

Once all the primes in the list are iterated and prime_tab is set accordingly, the newly discovered primes must be recorded. This occurs in the following loop:

register integer prime = ZERO;
for(register integer i = ZERO; i < UNIT; ++i) {
       	if(prime_tab[i] == false) {
               	prime = i + lstart;
               	prime_list->push_back(prime);
               	if(prime > max) {
                       	return;
               	}
       	}
}

Where the prime_tab buffer is iterated and all slots marked false (meaning they are prime) are used to compute the actual primes value which is the actual slot number + the logical start (i + lstart). This final computation actually generates the newly discovered prime. In the example above, slot 0 (prime_tab[0]) was marked composite and slot 2 was marked composite. However, slot 1 was marked prime. Therefore 1 + 4 == 5 which is prime. Note that for slot 3, it would have been marked prime as 3 + 4 == 7. At the end, if we have not discovered a prime equal to or larger than the user supplied max integer, we continue and repeat the process of prime discovery.

Using the Code

Using the code requires calling the C++ function indicating the maximum valued prime to generate and a pointer to a linked list which will contain the primes after the sieve is finished. An example usage is:

#include "bsieve.h"
int
main()
{
    list_int_t primes;
    if(argc != 2) {
        std::cerr << "syntax:: bsieve max_sized_prime" << std::endl;
        return 1;
    } 
    bsieve(strtoul(argv[1],NULL,10), &primes);   
    for(list_int_t::iterator i = primes.begin(); i != primes.end(); ++i) {
        std::cout << (*i) << std::endl;
    }
    return 0;
}

The code has been compiled and tested using g++ on Linux and MinGW on Windows.

Points of Interest

The largest number of primes that the algorithm can generate is dependent on the length of the linked list, so there is a limitation but hopefully this is better than simply allocating a large fixed size buffer. Further, any array allocation necessarily has half of its elements as composite as these are all even numbers.

The code as currently implemented uses an internal buffer called prime_tab which is allocated as static in thread local storage. It is possible to redesign the algorihtm to have this pointer supplied externally via the function arguments. Further, thread local storage is declared using compiler specific syntax.

At the point of insertion into the prime number list (prime_list), we can provide a callback to immediately try if a given integer is prime by testing if a newly discovered prime factors the integer. This can also be used to build a factoring algorithm for small integers. Future direction of the code should look into developing an algorithm that works in parallel using threads to generate primes and see if this technique can be utilized.

The code was tested on Linux and MinGW for Windows. The following bash loop was used to verify that all numbers generated are prime:

MAX_PRIME=1000;
for f in `./bsieve ${MAX_PRIME}`
do 
    factor $f  # All numbers should be prime.
done

References

• [1] Crandall, R. and Pomerance, C. (2005). Prime Numbers. A Computational Perspective. Second Edition. NY, NY: Springer