Eratosthenes 位阵列筛

Question

我正在尝试使用带有位数组的 Eratosthenes 筛来查找素数，但我使用的是 unsigned int 数组。我需要能够生成多达 2,147,483,647 个素数。我的代码可以工作并且可以生成大约 10,000,000，但是当我增加数组的大小以容纳更大的数字时，它会失败。有人可以指导我如何将位向量与 c（不是 c++）一起使用。 谢谢

这是我的代码：

#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;
}

Answer 1

使用位访问整数的示例
注意GetBit() 和SetBit()。
优化编译器将使/ 和% 快​​速使用2 的幂。

#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);
}

重写可以采用通常的建议不保存偶数，将 2 视为特殊情况。这有帮助，但我认为这是一个练习。通过使用 1 个字节将每 30 个倍数编码为 30 之后，我可以节省大约 4 倍，因为每 30 个整数最多有 8 个素数。存在其他方案。

Answer 2

我绝对不会使用位数组，而是使用原生 int 数组（取决于架构的 64 位或 32 位）并包装一个函数以将正常数字重新映射到正确的位置并按位 @987654321 @ 和&。

还要考虑省略偶数，它们几乎都不是素数。因此，您可以将前 128 个数字存储在第一个 64 位数字中，接下来的 128 个数字存储在第二个等中。

听起来有点复杂，但做起来还是挺有趣的！

Project Euler 似乎已经产生了一些非常好的解决方案。

好处是：筛分你不需要重新计算奇偶转移，但你可以取消设置每第三位筛分 3，每 5 位筛分 5 等等。

如果您想要快速的 java 解决方案作为详细参考，请来聊天。

EDIT4：更正了工作代码，但速度很慢。备忘：记得使用 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);
        }
        */
    }
}