Search code examples
cbitsieve-of-eratosthenesbitarray

Sieve of Eratosthenes bit array


I am trying to find primes using Sieve of Eratosthenes with bit arrays, but I am using unsigned int array. I need to be able to generate upto 2,147,483,647 primes. My code works and can generate around 10,000,000 but when I increase the size of my array to accommodate larger numbers it fails. Can someone guide me on how to use bit vectors with c (not c++). Thanks

Here's my code:

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

#define MAXBYTES 2000000
#define MAX 50000000
#define BITSIZE 32

void ClearBit(unsigned int [], unsigned int);
void SetBit(unsigned int [], unsigned int);
int  BitVal(unsigned int [], unsigned int);
void PrintBitStream(unsigned int [], unsigned long);
void PrintBitStreamData(unsigned int[], unsigned long);
int  Sieve(unsigned int[], unsigned int, unsigned int);

int main(int argc, char ** argv) {
    unsigned int maxsize = MAX;
    unsigned int i;
    //Set Bit Array
    unsigned int BitArray[MAXBYTES] = {0};
    SetBit(BitArray, 0);
    SetBit(BitArray, 1);
    i = 2;
    for (;i < maxsize;i++){

        if(Sieve(BitArray, i, maxsize)==0)
            break;
    }
    PrintBitStreamData(BitArray, maxsize-1);
    return EXIT_SUCCESS;
}

void PrintBitStreamData(unsigned int BitArray[], unsigned long maxsize) {
    unsigned int i;
    for (i = 0; i < maxsize; i++)
        if (!BitVal(BitArray, i))
            printf("%ld ", i);
    printf("\n");
}

void PrintBitStream(unsigned int BitArray[], unsigned long maxsize) {
    unsigned int i;
    for (i = 2; i < maxsize; i+=2)
        printf("%d", BitVal(BitArray, i));
    printf("\n");
}

void SetBit(unsigned int BitArray[], unsigned int pos) {
    BitArray[pos / BITSIZE] |= 1 << (pos % BITSIZE);
}

void ClearBit(unsigned int BitArray[], unsigned int pos) {
    BitArray[pos / BITSIZE] &= ~(1 << (pos % BITSIZE));
}

int BitVal(unsigned int BitArray[], unsigned int pos) {
    return ((BitArray[pos / BITSIZE] & (1 << (pos % BITSIZE))) != 0);
}

int Sieve(unsigned int BitArray[], unsigned int p, unsigned int maxsize) {
    unsigned int i;
    unsigned int j;
    j = 0;
    for (i = 2 * p; i < maxsize; i += p) {
        SetBit(BitArray, i);
        j++;
    }
    return j;
}

Solution

  • I definitely would NOT use a bit array, but an array of the native int (64-bit or 32-bit dependign on architecture) and wrapping around a function to remap normal numbers to the right place and bit with bitwise | and &.

    Also consider leaving out the even numbers, almost none of them are prime. So you could store the first 128 numbers in the first 64-bit number, the next 128 in the second etc.

    It sounds a bit complicated, but is a bit fun to get it work!

    Project Euler seems to have produced some really nice solutions.

    The good thing is: to sieve you don't need to recalculate the even-odd-transfer, but you can unset every third bit for sieving 3, every 5th bit for sieving 5 etc.

    Come into chat if you want a quick java solution as a detailed reference.

    EDIT4: corrected Working code, but yet slow. Memo: Remember using calloc!

    #include <stdio.h>
    #include <math.h>
    #include <stdlib.h>
    #include <limits.h>
    #include <time.h>
    
    typedef unsigned long long number;
    
    number lookFor = 2147483648ULL;
    number max = 2147483648ULL*10ULL; // hopefully more then every 10th uneven number is prime
    
    unsigned long * isComposite;
    
    number bitslong = CHAR_BIT*sizeof(long);
    
    time_t rawtime;
    struct tm * timeinfo;
    char buffer[80];
    
    // is not called during sieve, only once per sieving prime 
    // and needed for reading if a number is prime
    inline long getParts(number pos, number *slot, unsigned int *bit){
        *slot = pos / bitslong;
        *bit = (unsigned int)(pos % bitslong);
    }
    
    int isPrime(number n){
        if(n == 1){
            return 0;
        }
    
        if(n < 4){
            return 1;
        }
    
        if((n%2) == 0){
            return 0;
        }
    
        number slot=0;
        unsigned int bit=0;
        number pos = (number)(n-3)/2;
        getParts(pos, &slot, &bit);
        // printf("? n=%lld  pos = %lld slot = %lld bit = %lu ret %d \n", n, pos, slot, bit, !(isComposite[slot] & (1<<bit)));
        return !(isComposite[slot] & (1UL<<bit));
    }
    
    // start sieving at offset (internal position) offset with width step 
    int doSieve(number offset, number step){
        rawtime = time(0);
        time (&rawtime);
        timeinfo = localtime (&rawtime);
    
        strftime(buffer, 80, "%Y-%m-%d %H:%I:%S", timeinfo);
        unsigned int bit=0;
        number slot=0;
        getParts(offset, &slot, &bit);
        printf("doSieve %s  %lld %lld  %lu \n", buffer, offset, step, isComposite[slot]);
    
        while(slot < max/bitslong){
            slot += (step + bit)/bitslong;
            bit = (step + bit) % bitslong;
            isComposite[slot] |= (1UL << bit);
        } 
        return 1;
    }
    
    int sieve(){
        number spot;
        spot=1;
        number pos;
        pos = 0;
        while(spot < 1 + sqrt((float)max)){
            spot+=2;
            if(! isPrime(spot)){
                pos++;
                continue;
            }
            doSieve(pos, spot);
            pos++;
        }
    }
    
    void main(int argc, char *argv[]){
        if(argc > 1){
            char *tp = malloc(sizeof(char*));
            max = strtol(argv[1], &tp, 10);
        }
        printf("max %lld , sq %ld, malloc: %lld\n", max, (long)(1 + sqrt((float)max)), 1+max/bitslong);
        isComposite = calloc((2+max/bitslong), sizeof(unsigned long)) ;
        if(! isComposite){
            printf("no mem\n");
            exit(5);
        }
        sieve();
        number i;
        number found = 0;
        for(i = 1; i<max && found < lookFor; i++){
            if(isPrime(i)){
                found++;
                // printf(" %30lld %30lld \n", found, i);
                if(found % 10000 == 0 ){
                    printf("%30lld %30lld \n", found, i);
                }
            }
            /*
            if(i % 1000 == 17){
                printf("%5lld %5lld \n", i, found);
            }
            */
        }
    }