Counting number of primes within a given range of long int using C

Question

Well, there are lots of such questions available in SO as well as other forums. However, none of these helped.

I wrote a program in "C" to find number of primes within a range. The range i in long int. I am using Sieve of Eratosthenes" algorithm. I am using an array of long ints to store all the numbers from 1 till the limit. I could not think of a better approach to achieve without using an array. The code works fine, till 10000000. But after that, it runs out of memory and exits. Below is my code.

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
typedef unsigned long uint_32;

int main() {

    uint_32 i, N, *list, cross=0, j=4, k, primes_cnt = 0;
    clock_t start, end;
    double exec_time;

    system("cls");
    printf("Enter N\n");
    scanf("%lu", &N);
    list = (uint_32 *) malloc( (N+1) * sizeof(uint_32));
    start = clock();
    for(i=0; i<=N+1; i++) {
        list[i] = i;
    }

    for(i=0; cross<=N/2; i++) { 
        if(i == 0)
            cross = 2;
        else if(i == 1)
            cross = 3;
        else {

            for(j=cross+1; j<=N; j++) {
                if(list[j] != 0){

                    cross = list[j];
                    break;
                }
            }
        }

        for(k=cross*2; k<=N; k+=cross) {
            if(k <= N)
                list[k] = 0;
        }
    }

    for(i=2; i<=N; i++) {

        if(list[i] == 0)
            continue;
        else
            primes_cnt++;
    }

    printf("%lu", primes_cnt);
    end = clock();
    exec_time = (double) (end-start);
    printf("\n%f", exec_time);
    return 0;
}

I am stuck and can't think of a better way to achieve this. Any help will be hugely appreciated. Thanks.

Edit:

My aim is to generate and print all prime numbers below the range. As printing consumed a lot of time, I thought of getting the first step right.

Answer 1

I could see that the approach you are using is the basic implementation of Eratosthenes, that first stick out all the 2's multiple and then 3's multiple and so on.

But I have a better solution to the question. Actually, there is question on spoj PRINT. Please go through it and do check the constraints it follows. Below is my code snippet for this problem:

#include<stdio.h>
#include<math.h>
#include<cstdlib>

int num[46500] = {0},prime[5000],prime_index = -1;

int main() {
    /* First, calculate the prime up-to the sqrt(N) (preferably greater than, but near to 
       sqrt(N) */   
    prime[++prime_index] = 2;    int i,j,k;

    for(i=3;    i<216;     i += 2) {
        if(num[i] == 0) {
            prime[++prime_index] = i;
            for(j = i*i, k = 2*i;    j<=46500;   j += k) {
                num[j] = 1;
            }
        }
    }

    for(;   i<=46500;   i+= 2) {
        if(num[i] == 0) {
            prime[++prime_index] = i;
        }
    }

    int t;            // Stands for number of test cases  

    scanf("%i",&t);
    while(t--) {

        bool arr[1000005] = {0};   int m,n,j,k;

        scanf("%i%i",&m,&n);

        if(m == 1)
            m++;

        if(m == 2 && m <= n) {
            printf("2\n");
        }

        int sqt = sqrt(n) + 1;

        for(i=0;    i<=prime_index; i++) {
            if(prime[i] > sqt) {
                sqt = i;
                break;
            }
        }

        for(;   m<=n && m <= prime[prime_index];   m++) {
            if(m&1 && num[m] == 0) {
                printf("%i\n",m);
            }
        }

        if(m%2 == 0) {
            m++;
        }

        for(i=1;    i<=sqt;     i++) {
            j = (m%prime[i]) ? (m + prime[i] - m%prime[i]) : (m);
            for(k=j;    k<=n;   k += prime[i]) {
                arr[k-m] = 1;
            }
        }

        for(i=0;    i<=n-m; i += 2) {
            if(!arr[i]) {
                printf("%i\n",m+i);
            }
        }
        printf("\n");
    }
    return 0;
}

I hope you got the point:

And, as you mentioned that your program is working fine up-to 10^7 but above it fails, it must be because you must be running out of the memory.

NOTE: I'm sharing my code only for knowledge purpose. Please, don't copy and paste it, until you get the point.

Answer 2

You can tweak the algorithm a bit to calculate the prime numbers in chunks.

Load a part of the array (as much as fits the memory), and in addition hold a list of all known prime numbers.
Whenever you load a chunk, first go through the already known prime numbers, and similar to the regular sieve, set all non primes as such.
Then, go over the array again, mark whatever you can, and add to the list the new prime numbers found.

When done, you'll have a list containing all your prime numbers.

Answer 3

There are other algorithm that does not require you to generate prime number up to N to count number of prime below N. The easiest algorithm to implement is Legendre Prime Counting. The algorithm requires you to generate only sqrt(N) prime to determine the number of prime below N.

The idea behind the algorithm is that

pi(n) = phi(n, sqrt(n)) + pi(sqrt(n)) - 1

where
   pi(n) = number of prime below N
   phi(n, m) = number of number below N that is not divisible by any prime below m.

That's mean phi(n, sqrt(n)) = number of prime between sqrt(n) to n. For how to calculate the phi, you can go to the following link (Feasible implementation of a Prime Counting Function)

The reason why it is more efficient is because it is easiest to compute phi(n, m) than to compute pi(n). Let say that I want to compute phi(100, 3) means that how many number below or equal to 100 that does not divisible by 2 and 3. You can do as following. phi(100, 3) = 100 - 100/2 - 100/3 + 100/6.

Answer 4

Your code uses about 32 times as much memory as it needs. Note that since you initialized list[i] = i the assignment cross = list[j] can be replaced with cross = j, making it possible to replace list with a bit vector.

However, this is not enough to bring the range to 2⁶⁴, because your implementation would require 2⁶¹ bytes (2 exbibytes) of memory, so you need to optimize some more.

The next thing to notice is that you do not need to go up to N/2 when "crossing" the numbers: √N is sufficient (you should be able to prove this by thinking about the result of dividing a composite number by its divisors above √N). This brings memory requirements within your reach, because your "crossing" primes would fit in about 4 GB of memory.

Once you have an array of crossing primes, you can build a partial sieve for any range without keeping in memory all ranges that precede it. This is called the Segmented sieve. You can find details on it, along with a simple implementation, on the page of primesieve generator. Another advantage of this approach is that you can parallelize it, bringing the time down even further.