Listing prime numbers using Sieve's method using bitmask

Question

I wrote the following code to list all the prime numbers upto 2 billion using Sieve's method. I used bitmasking for flagging purpose. While I am able to get the prime numbers correctly, a few primes in the beginning are missing every time. Please help me find the bug in the program.

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

#define MAX 2000000000

char* listPrimes(){
int block = sqrt(MAX);
char* mark = calloc((MAX/8),sizeof(char));
int i = 2;
int j;
char mask[8];
for(j=0;j<8;j++)
    mask[j] = 0;
mask[7] = 1;
mask[6] |= mask[7] << 1;
mask[5] |= mask[7] << 2;
mask[4] |= mask[7] << 3;
mask[3] |= mask[7] << 4;
mask[2] |= mask[7] << 5;
mask[1] |= mask[7] << 6;
mask[0] |= mask[7] << 7;

for(j=0;j<8;j++)
    printf("%d ",mask[j]);
mark[0] |= mask[0];
mark[0] |= mask[1];

while (i < block){

        for (j = 2; i*j <= block; j++)
                mark[(i*j) / 8] |= mask[((i*j) % 8 )];
        i++;
    }
printf("\n");
printf("The block size is\t:\t%d\n",block);


j = 2;
while(j<=block){
    if((mark[j / 8] & mask[j]) == 0 ){
        for(i = 2;i <= MAX; i++){
            if((i%j) == 0){
                mark[i / 8] |= mask[(i % 8)];
            }
        }
    }
while((mark[++j / 8] & mask[j % 8]) != 0);
}


for(j=0;j<=MAX;j++)
        if((mark[j / 8] & mask[(j % 8)]) == 0)
            printf("%d\n", ((8*(j / 8)) + (j % 8)));

return mark;
}   

int main(int argc,char* argv[]){

listPrimes();

return 0;
}

Answer 1

As ArunMK said, in the second while loop you mark the prime j itself as a multiple of j. And as Lee Meador said, you need to take the modulus of j modulo 8 for the mask index, otherwise you access out of bounds and invoke undefined behaviour.

A further point where you invoke undefined behaviour is

while((mark[++j / 8] & mask[j % 8]) != 0);

where you use and modify j without intervening sequence point. You can avoid that by writing

do {
    ++j;
}while((mark[j/8] & mask[j%8]) != 0);

or, if you insist on a while loop with empty body

while(++j, (mark[j/8] & mask[j%8]) != 0);

you can use the comma operator.

More undefined behaviour by accessing mark[MAX/8] which is not allocated in

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

and

for(j=0;j<=MAX;j++)

Also, if char is signed and eight bits wide,

mask[0] |= mask[7] << 7;

is implementation-defined (and may raise an implementation-defined signal) since the result of

mask[0] | (mask[7] << 7)

(the int 128) is not representable as a char.

But why are you dividing each number by all primes not exceeding the square root of the bound in the second while loop?

    for(i = 2;i <= MAX; i++){
        if((i%j) == 0){

That makes your algorithm not a Sieve of Eratosthenes, but a trial division.

Why don't you use the technique from the first while loop there too? (And then, why two loops at all?)

while (i <= block){
    if ((mark[i/8] & mask[i%8]) == 0) {
        for (j = 2; i*j < MAX; j++) {
            mark[(i*j) / 8] |= mask[((i*j) % 8 )];
        }
    }
    i++;
}

would not overflow (for the given value of MAX, if that is representable as an int), and produce the correct output orders of magnitude faster.

Answer 2

Change the middle loop to add the modulo:

j = 2;
while(j<=block){
    if((mark[j / 8] & mask[j % 8]) == 0 ){
        for(i = 2;i <= MAX; i++){
            if((i%j) == 0){
                mark[i / 8] |= mask[(i % 8)];
            }
        }
    }
}

Answer 3

In the second while loop you are looping through i from 2 onwards and you do an if (i%j == 0). This will be true for i when it is a prime as well. You need to check for (i != j). Also the modulo as reported above. Hence it becomes: if ((i%j == 0) { if (i!=j) mark[i/j] |= mask[i%j]; }