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

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.

• This one is better. – Ragib Zaman Feb 21 '12 at 10:38
• I would suggest you wheel factorization: en.wikipedia.org/wiki/Wheel_factorization – rubik Feb 21 '12 at 12:40
• Why does the title talk about random numbers? Am I missing something? – Peter Taylor Feb 21 '12 at 15:06
• If your program ever crashes then it definitely needs fixing. – Nate Eldredge Feb 21 '12 at 15:49
• random numbers? – Marty Feb 21 '12 at 23:26

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;   // primes to use in the sieve
int First;    // 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.

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:

• it's an assignment for class. i think that would be plagiarism :) – Marty Feb 22 '12 at 0:39

I came up with a different wheel sieve ( not coded ) and would be glad to send it to you. It's a bit lengthy to post here. my email is madjatgoldendotnet. (Sorry for all the sins I commited with this answer )

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".

• but i still need to iterate up to 2^32 to get rid of those multiples. – Marty Feb 21 '12 at 23:27
• @Frederick: There are $4,091,687,074$ composite numbers up to $2^{32}$. If removing those from a list is possible in some sense, then about $2.8$ times as many steps, as required by the sieve of Eratosthenes, should usually also be possible. – Henry Feb 22 '12 at 1:03
• I must have something wrong, it's been going for ten minutes with 8 threads going and it's only at multiples of 8537. guess I'll hafta head over to stackoverflow... – Marty Feb 22 '12 at 5:57