Question

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;
}
Was it helpful?

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);
        }
        */
    }
}

OTHER TIPS

Example using bit access into an integer
Note GetBit() and SetBit().
An optimizing compiler will make the / and % fast for powers of 2 are used.

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

#define ubitsize (sizeof(unsigned)*8)

unsigned GetBit(const unsigned *List, unsigned long index) {
  return !!(List[index / ubitsize] & (1u << (index % ubitsize)));
  }

void SetBit(unsigned *List, unsigned long index) {
  List[index / ubitsize] |= (1u << (index % ubitsize));
  }

void Sieve_of_Eratosthenes_via_bit_array(unsigned long MaxCandidatePrime) {
  unsigned long uByteSize = MaxCandidatePrime/ubitsize + 1;
  unsigned *List = calloc(uByteSize, sizeof *List);
  if (List == 0) return;

  unsigned long PrimeCount = 0;
  unsigned long Index = 0;
  for (Index = 2; Index <= MaxCandidatePrime; Index++) {
    // If found a prime ...
    if (GetBit(List, Index) == 0) {
      PrimeCount++;
      // let's see the progress
      if ((PrimeCount % (256LU*1024)) == 0) printf("%lu\n", Index);

      // Mark subsequent multiples as not--a-prime
      unsigned long Index2 = Index*2;
      while (Index2 <= MaxCandidatePrime) {
        SetBit(List, Index2);
        Index2 += Index;
      }
    }
  }
  printf("X %lu\n", Index);
  free(List);
}

void test(void) {
  Sieve_of_Eratosthenes_via_bit_array(200LU*1000*1000);
}

A re-write could employ usual suggestion to not save even numbers, treating 2 as a special case. This helps but I assume this is an exercise. I could save a factor of about 4 by saving using 1 byte to encode every 30th multiple as, after 30, for there are at most 8 primes every 30 integers. Other schemes exist.

Licensed under: CC-BY-SA with attribution
Not affiliated with StackOverflow
scroll top