Listing prime numbers up to 2 billion using sieve's method in C

Question

I tried listing prime numbers up to 2 billion, using Sieve Eratosthenes method. Here is what I used!

The problem I am facing is, I am not able to go beyond 10 million numbers. When I tried, it says 'Segmentation Fault'. I searched in the Internet to find the cause. Some sites say, it is the memory allocation limitation of the compiler itself. Some say, it is a hardware limitation. I am using a 64-bit processor with 4GB of RAM installed. Please suggest me a way to list them out.

#include <stdio.h>
#include <stdlib.h>

#define MAX 1000000
long int mark[MAX] = {0};

int isone(){
   long int i;
   long int cnt = 0;
   for(i = 0; i < MAX ; i++){
      if(mark[i] == 1)
         cnt++;
   }
   if(cnt == MAX)
      return 1;
   else
      return 0;
}



int main(int argc,char* argv[]){
   long int list[MAX];
   long int s = 2;
   long int i;
   for(i = 0 ; i < MAX ; i++){
      list[i] = s;
      s++;
   }
   s = 2;
   printf("\n");

   while(isone() == 0){
      for(i = 0 ; i < MAX ; i++){
         if((list[i] % s) == 0)
            mark[i] = 1;
      }
      printf(" %lu ",s);
      while(mark[++s - 2] != 0);
   }

   return 1;
}

Answer 1

long int mark[1000000] does stack allocation, which is not what you want for that amount of memory. try long int *mark = malloc(sizeof(long int) * 1000000) instead to allocate heap memory. This will get you beyond ~1Mil of array elements.

remember to free the memory, if you don't use it anymore. if yon don't know malloc or free, read the manpages (manuals) for the functions, available via man 3 malloc and man 3 free on any linux terminal. (alternatively you could just google man malloc)

EDIT: make that calloc(1000000, sizeof(long int)) to have a 0-initialized array, which is probably better.

Additionally, you can use every element of your array as a bitmask, to be able to store one mark per bit, and not per sizeof(long int) bytes. I'd recommend using a fixed-width integer type, like uint32_t for the array elements and then setting the (n % 32)'th bit in the (n / 32)'th element of the array to 1 instead of just setting the nth element to 1.

you can set the nth bit of an integer i by using:

uint32_t i = 0;
i |= ((uint32_t) 1) << n

assuming you start counting at 0.

that makes your set operation on the uint32_t bitmask array for a number n:

mask[n / 32] |= ((uint32_t)1) << (n % 32)

that saves you >99% of memory for 32bit types. Have fun :D

Another, more advanced approach to use here is prime wheel factorization, which basically means that you declare 2,3 and 5 (and possibly even more) as prime beforehand, and use only numbers that are not divisible by one of these in your mask array. But that's a really advanced concept.

However, I have written a primesieve wich wheel factorization for 2 and 3 in C in about ~15 lines of code (also for projecteuler) so it is possible to implement this stuff efficiently ;)

Answer 2

The most immediate improvement is to switch to bits representing the odd numbers. Thus to cover the M=2 billion numbers, or 1 billion odds, you need 1000/8 = 125 million bytes =~ 120 MB of memory (allocate them on heap, still, with the calloc function).

The bit at position i will represent the number 2*i+1. Thus when marking the multiples of a prime p, i.e. p^2, p^2+2p, ..., M, we have p^2=(2i+1)^2=4i^2+4i+1 represented by a bit at the position j=(p^2-1)/2=2i(i+1), and next multiples of p above it at position increments of p=2i+1,

for( i=1; ; ++i )
  if( bit_not_set(i) )
  {
      p=i+i+1;
      k=(p-1)*(i+1);
      if( k > 1000000000) break;
      for( ; k<1000000000; k+=p)
         set_bit(k); // mark as composite
  }
// all bits i>0 where bit_not_set(i) holds, 
// represent prime numbers 2i+1

Next step is to switch to working in smaller segments that will fit in your cache size. This should speed things up. You will only need to reserve memory region for primes under the square root of 2 billion in value, i.e. 44721.

First, sieve this smaller region to find the primes there; then write these primes into a separate int array; then use this array of primes to sieve each segment, possibly printing the found primes to stdout or whatever.