Is it a bad idea to use a Sieve of Eratosthenes to find all primes up to very large N?

Question

I need to write a program in C++ that finds all primes up to 2^32. I used a Sieve of Eratosthenes with multiple threads, but it only worked well up to about 10 million. After that it just takes too long. Sometimes it would even crash. I'm wondering if it's my code that needs fixing, or if it's just a bad idea to use Sieve of Eratosthenes with that big a number? What are better prime sieves for large numbers? I've read about some other ones, but can't find anywhere that says which are good for larger ranges.

Answer 1

The sieve is very sensible to memory use, and it's speed depends a lot on the implementation, if you pick the "classical":

bool C[1u<<32];  // this is too big
void sieve(unsigned Lim){
    for( unsigned p = 2; p*p < Lim; ++p)
        if( !C[p] )
            for( unsigned n = p*p; n < Lim; n+=p )
                C[n] = true;
}

then it does not even start running in most system due to lack of space to hold C. even if you compress the integer info in bits (you can keep up to 30 integers / byte using mod 30 and discarding multiples of 2,3 and 5) the resulting array is very big and inefficient.

You can multiply the efficiency by a couple of orders of magnitude processing small blocks (say from 60k to 1MB depending in your computer's cache size), the idea is to pre-compute the primes used in the sieve and to keep for each one the next multiple, as in the following code:

const int BLOCKSIZE =1<<16;
bool C[BLOCKSIZE];
int Primes[10000];   // primes to use in the sieve
int First[10000];    // next integer divisible by p (delta form block start)
int NP = 0;

void sieve(unsigned Lim){
    for( unsigned p = 2; (long long)p*p <= Lim; ++p)
        if( !C[p] ){
            Primes[NP] = p;
            First[NP] = p*p;
            NP++;
            for( unsigned n = p*p; (long long)n*n <= Lim; n+=p )
                C[n] = 1;
        }

    for( unsigned B = 0; B < Lim; B+=BLOCKSIZE ){
        for(int j = 0; j < BLOCKSIZE; j++)
            C[j] = 0;
        for(int k = 0; k < NP; ++k){
            int p = Primes[k];
            int t = First[k];
            for( ; t < BLOCKSIZE; t+=p ){
                C[t] = 1;
            }
            First[k] = t-BLOCKSIZE;
        }
        // Block is finished, primes are the integers B+m with C[m]==0.
        ...   // process the primes here
    }
}

There are a few details missing, and quite a few other improvements you can make in this code, compress 30 integers in 8 bits, use wheels, ... .

There are sieves that are, at least theoretically, faster than the sieve of Eratosthenes, namely the sieve of Atkin, but it is very difficult to implement and probably it is slower for this order of magnitude.

Answer 2

The Sieve of Eratosthenes is very good in that range, so it's probably your implementation. But why roll your own? There are excellent implementations out there:

• http://sites.google.com/site/bbuhrow/
• http://code.google.com/p/primesieve/

Answer 3

If you are looking for primes up to $2^{32}$ then you only need to sieve out multiples of primes up to $2^{16} = 65536$.

All the remaining numbers up to $2^{32}$ (apart from 1) will be prime.

So you do not need to look "up to about 10 million".