Prime-sieve optimized for memory

Question

I am looking for a prime-sieve implementation which is efficient in terms of memory consumption.

Of course, the primality-test itself should execute at a constant and minimal number of operations.

I have implemented a sieve that indicates primality only for numbers that are adjacent to multiples of 6.

For any other number, either it is 2 or 3 (therefore prime), or it is a multiple of 2 or 3 (therefore non-prime).

So this is what I came up with, and I've been wondering if there is anything better under those requirements:

Interface:

#include <limits.h>

// Defined by the user (must be less than 'UINT_MAX')
#define RANGE 4000000000

// The actual length required for the prime-sieve array
#define ARR_LEN (((RANGE-1)/(3*CHAR_BIT)+1))

// Assumes that all entries in 'sieve' are initialized to zero
void Init(char sieve[ARR_LEN]);

// Assumes that 'Init(sieve)' has been called and that '1 < n < RANGE'
int IsPrime(char sieve[ARR_LEN],unsigned int n);

#if RANGE >= UINT_MAX
    #error RANGE exceeds the limit
#endif

Implementation:

#include <math.h>

#define GET_BIT(sieve,n) ((sieve[(n)/(3*CHAR_BIT)]>>((n)%(3*CHAR_BIT)/3))&1)
#define SET_BIT(sieve,n) sieve[(n)/(3*CHAR_BIT)] |= 1<<((n)%(3*CHAR_BIT)/3)

static void InitOne(char sieve[ARR_LEN],int d)
{
    unsigned int i,j;
    unsigned int root = (unsigned int)sqrt((double)RANGE);

    for (i=6+d; i<=root; i+=6)
    {
        if (GET_BIT(sieve,i) == 0)
        {
            for (j=6*i; j<RANGE; j+=6*i)
            {
                SET_BIT(sieve,j-i);
                SET_BIT(sieve,j+i);
            }
        }
    }
}

void Init(char sieve[ARR_LEN])
{
    InitOne(sieve,-1);
    InitOne(sieve,+1);
}

int IsPrime(char sieve[ARR_LEN],unsigned int n)
{
    return n == 2 || n == 3 || (n%2 != 0 && n%3 != 0 && GET_BIT(sieve,n) == 0);
}

Answer 1

You've correctly deduced that you can exploit the fact that there are only two number that are relatively prime to 6, i.e. 1 and 5 (aka +1 and -1). Using that fact and storing the sieve as bits instead of bytes, you reduce the memory requirement by a factor of 24.

To save a little more memory, you can go to the next level and note that there are only 8 numbers (modulo 30) that are relatively prime to 30. They are 1, 7, 11, 13, 17, 19, 23, 29. Using that fact and storing as bits, the memory is reduced by a factor of 30.

Implementation note: each byte in the sieve represents the 8 numbers relatively prime to some multiple of 30. For example, the bits contained in sieve[3] represent the numbers 91, 97, 101, ...

Answer 2

For numbers less than the max bit size you can use a single "primes" mask with a 1 in each prime position.

if (n < 64) return !!(primes & 1ull<<n);

To filter out multiples of 2, 3 & 5 for higher values of n you can use n%30 along with a notprime mask (note that the product of the first 3 primes is a constant of 30, so most compilers will optimize it to multiply and shift)

else if (notprime & 1ull << n%30) return 0;

For the remaining values that may be prime, you can do a manual check.

else return isprime(n);

Since each block of 30 numbers contains 8 possible primes IIRC, the isprime function can be vectorized by simply adding 30 to each successive block or you can trade size for speed and store precalculated multiply and shift values in a constant vector array to avoid expensive divisions. Alternatively you can use 8bit masks to store which of the maybe-primes are prime as you go, so that you only need to compute it once.