Sieve of Eratosthenes - bitwise optimization problem

Question

I guess every one of you has encountered the optimization code of Eratosthenes sieve with bitwise operations. I am trying to wrap my head around it and I have a question on one of the operations in this implementation. Here is the code from GeeksforGeeks:

bool ifnotPrime(int prime[], int x) { 
    // checking whether the value of element 
    // is set or not. Using prime[x/64], we find 
    // the slot in prime array. To find the bit 
    // number, we divide x by 2 and take its mod 
    // with 32. 
    return (prime[x / 64] & (1 << ((x >> 1) & 31))); 
} 

// Marks x composite in prime[] 
bool makeComposite(int prime[], int x) { 
    // Set a bit corresponding to given element. 
    // Using prime[x/64], we find the slot in prime  
    // array. To find the bit number, we divide x 
    // by 2 and take its mod with 32. 
    prime[x / 64] |= (1 << ((x >> 1) & 31)); 
} 

// Prints all prime numbers smaller than n. 
void bitWiseSieve(int n) { 
    // Assuming that n takes 32 bits, we reduce 
    // size to n/64 from n/2. 
    int prime[n / 64]; 

    // Initializing values to 0 . 
    memset(prime, 0, sizeof(prime)); 

    // 2 is the only even prime so we can ignore that 
    // loop starts from 3 as we have used in sieve of 
    // Eratosthenes . 
    for (int i = 3; i * i <= n; i += 2) { 

        // If i is prime, mark all its multiples as 
        // composite 
        if (!ifnotPrime(prime, i)) 
            for (int j = i * i, k = i << 1; j < n; j += k) 
                makeComposite(prime, j); 
    } 

    // writing 2 separately 
    printf("2 "); 

    // Printing other primes 
    for (int i = 3; i <= n; i += 2) 
        if (!ifnotPrime(prime, i)) 
            printf("%d ", i); 
} 

// Driver code 
int main() { 
    int n = 30; 
    bitWiseSieve(n); 
    return 0; 
} 

So my questions are:

1. what is the meaning of (prime[x/64] & (1 << ((x >> 1) & 31)) and more specifically (1 << ((x >> 1) & 31));
2. in prime[x/64] why do we divide by 64 and not with 32, when we work with 32 bit integer;
3. Is int prime[n/64] correct if n < 64?

Answer 1

There are multiple problems in the code:

• makeComposite() should have return type void.
• (1 << ((x >> 1) & 31)) has undefined behavior if x == 63 because 1 << 31 overflows the range of type int. You should use 1U or preferably 1UL to ensure 32 bits.
• it would be simpler to shift the table entry and mask the last bit: return (prime[x / 64] >> ((x >> 1) & 31)) & 1;
• instead of assuming that type int has 32 bits, you should use uint32_t as the type for the bit array.
• the array int prime[n / 64]; is too short if n is not a multiple of 64. Use uint32_t prime[n / 64 + 1]; instead. This is a problem for your example as n = 30, so the array is created with a length of 0.
• ifnotPrime(n) only returns a valid result for odd values. It might be better and more readable to change this function and name it isOddPrime().

Regarding your questions:

what is the meaning of (prime[x/64] & (1 << ((x >> 1) & 31)) and more specifically (1 << ((x >> 1) & 31))?

x is first divided by 2 (shift right one position) because only odd numbers have a bit in the array, then the result is masked with 31 to keep the 5 low order bits as the bit number in the word. For any unsigned value x, x & 31 is equivalent to x % 32. x / 64 is the word number where to test this bit.

As mentioned above, 1 is an int and thus should not be shifted left by 31 positions. Using 1UL ensures that the type has at least 32 bits and can be shifted by 31 positions.

in prime[x/64] why do we divide by 64 and not with 32, when we work with 32-bit integer?

The bits in the array correspond to odd numbers, so a 32-bit word contains primeness information for 64 numbers: 32 even numbers that are known to be composite and 32 odd numbers for which the corresponding bit it set if the number is composite.

Is int prime[n/64] correct if n < 64?

Not it is not, it is incorrect if n is not a multiple of 64: the size expression should be (n + 63) / 64, or better int prime[n/64 + 1]

---

Here is a modified version, where you can pass a command line argument:

#include <stdio.h>
#include <stdint.h>
#include <stdlib.h>
#include <stdbool.h>
#include <string.h>

bool isOddPrime(const uint32_t prime[], unsigned x) {
    // checking whether the value of element
    // is set or not. Using prime[x/64], we find
    // the slot in prime array. To find the bit
    // number, we divide x by 2 and take its mod
    // with 32.
    return 1 ^ ((prime[x / 64] >> ((x >> 1) & 31)) & 1);
}

// Marks x composite in prime[]
void makeComposite(uint32_t prime[], unsigned x) {
    // Set a bit corresponding to given element.
    // Using prime[x/64], we find the slot in prime
    // array. To find the bit number, we divide x
    // by 2 and take its mod with 32.
    prime[x / 64] |= (1UL << ((x >> 1) & 31));
}

// Prints all prime numbers smaller than n.
void bitWiseSieve(unsigned n) {
    // Assuming that n takes 32 bits, we reduce
    // size to n/64 from n/2.
    uint32_t prime[n / 64 + 1];

    // Initializing values to 0 .
    memset(prime, 0, sizeof(prime));

    // 2 is the only even prime so we can ignore that
    // loop starts from 3 as we have used in sieve of
    // Eratosthenes .
    for (unsigned i = 3; i * i <= n; i += 2) {
        // If i is prime, mark all its multiples as composite
        if (isOddPrime(prime, i)) {
            for (unsigned j = i * i, k = i << 1; j < n; j += k)
                makeComposite(prime, j);
        }
    }

    // writing 2 separately
    if (n >= 2)
        printf("2\n");

    // Printing other primes
    for (unsigned i = 3; i <= n; i += 2) {
        if (isOddPrime(prime, i))
            printf("%u\n", i);
    }
}

// Driver code
int main(int argc, char *argv[]) {
    unsigned n = argc > 1 ? strtol(argv[1], NULL, 0) : 1000;
    bitWiseSieve(n);
    return 0;
}