arrays - 优化 Eratosthenes 筛法，向 Bool 数组添加轮子

Question

我知道有几篇关于如何实现轮子的帖子，但我真的很想看看如何用我目前的筛子方法来实现一个轮子。

我在 C 中创建了自己的位数组，实现如下:

#define setBit1(Array, i) (Array[i/INT_BITS] |= (1 << i % INT_BITS))
#define getBit(Array, i) ((Array[i/INT_BITS] & (1 << i % INT_BITS)) ? 1 : 0)
#define setBit0(Array, i) (Array[i/INT_BITS] &= ~(1 << i % INT_BITS))

int * createBitArray(unsigned long long size) {

    // find bytes required, round down to nearest whole byte
    unsigned long long bytesRequired = size / BITS_PERBYTE;

    // round up to highest multiple of 4 or 8 bytes (one int) 
    bytesRequired = (sizeof(int) * (bytesRequired / sizeof(int) + 
        ((size % BITS_PERBYTE * sizeof(int) == 0) ? 0 : 1)));

    // allocate array of "bits", round number of ints required up 
    return (int *)malloc((bytesRequired)); 

}

我在 C 语言中使用 clock() 进行了一些测试，我发现对于超过 1,000,000 的大型数组，位数组，即使进行位操作，也至少比它快 200%一个整数数组。它还使用了 1/32 的内存。

#define indexToNum(n) (2*n + 1) 
#define numToIndex(n) ((n - 1) / 2)
typedef unsigned long long LONG; 

// populates prime array through Sieve of Eratosthenes, taking custom 
// odd keyed bit array, and the raw array length, as arguments
void findPrimes(int * primes, LONG arrLength) {

    long sqrtArrLength = (long)((sqrt((2 * arrLength) + 1) - 1) / 2); 
    long maxMult = 0; 
    long integerFromIndex = 0; 

    for (int i = 1; i <= sqrtArrLength; i++) {

         if (!getBit(primes, i)) {
             integerFromIndex = indexToNum(i); 
             maxMult = (indexToNum(arrLength)) / integerFromIndex; 

             for (int j = integerFromIndex; j <= maxMult; j+= 2) {
                 setBit1(primes, numToIndex((integerFromIndex*j)));     
            }

        }   
    }
}

我一直在用索引 i 填充位数组，代表通过 (2i + 1) 获得的数字。这样做的好处是减少了迭代偶数所花费的时间，并再次将数组的必要内存减少了一半。 2 被手动添加到素数之后。这是在索引和数字之间转换所花费时间的结果，但根据我的测试，对于超过 1,000 个素数，这种方法更快。

我对如何进一步优化感到困惑；我已经减小了数组大小，我只测试 sqrt(n)，我开始从 p * p 向上“筛选”素数，我已经消除了所有偶数，它仍然需要我大约 60 秒C 中的前 100,000,000 个素数。

据我所知，“wheel”方法要求将数字的实际整数存储在索引中。我真的坚持用我当前的位数组来实现它。

Answer 1

当我修复您的实现并在我的 Macbook Pro 上运行它时，标记所有复合 Material 需要 17 秒 <= 2^31，这非常快。

不过，您还可以尝试其他一些事情。使用轮子可能会减少一半的时间。

如果仔细实现，欧拉筛法是线性时间，但它需要一个 int 数组而不是位数组。

阿特金筛需要线性时间并且非常实用:https://en.wikipedia.org/wiki/Sieve_of_Atkin

最后是我自己的(这只是意味着我没有在其他任何地方看到它，但我也没有真正看过) super 有趣的筛子也需要线性时间并在 6.5 秒内找到所有 <= 2^31 的素数。感谢您给我发帖的借口:

#include <stdio.h>
#include <time.h>
#include <stdlib.h>
#include <memory.h>
#include <math.h>

#define SETBIT(mem,num) ( ((uint8_t *)mem)[(num)>>4] |= ((uint8_t)1)<<(((num)>>1)&7) )

int main(int argc, char *argv[])
{
    //we'll find all primes <= this
    uint32_t MAXTEST = (1L<<31)-1;
    //which will be less than this many
    size_t MAXPRIMES = 110000000;

    //We'll find this many primes with the sieve of Eratosthenes.
    //After that, we'll switch to a linear time algorithm
    size_t NUMPREMARK = 48;

    //Allocate a bit array for all odd numbers <= MAXTEST
    size_t NBYTES = (MAXTEST>>4)+1;
    uint8_t *mem = malloc(NBYTES);
    memset(mem, 0, NBYTES);

    //We'll store the primes in here
    unsigned *primes = (unsigned *)malloc(MAXPRIMES*sizeof(unsigned));
    size_t nprimes = 0;


    clock_t start_t = clock();

    //first premark with sieve or Eratosthenes
    primes[nprimes++]=2;
    for (uint32_t test=3; nprimes<NUMPREMARK; test+=2)
    {
        if ( mem[test>>4] & ((uint8_t)1<<((test>>1)&7)) )
        {
            continue;
        }
        primes[nprimes++]=test;
        uint32_t inc=test<<1;
        for(uint32_t prod=test*test; prod<=MAXTEST; prod+=inc)
        {
            SETBIT(mem,prod);
        }
    }


    //Iterate through all products of the remaining primes and mark them off, 
    //in linear time. Products are visited in lexical order of their  
    //prime factorizations, with factors in descending order.

    //stacks containing the current prime indexes and partial products for 
    //prefixes of the current product
    size_t stksize=0;
    size_t indexes[64];
    uint32_t products[64];

    for (uint32_t test=primes[NUMPREMARK-1]+2; test<=MAXTEST; test+=2)
    {
        if ( mem[test>>4] & ((uint8_t)1<<((test>>1)&7)) )
        {
            continue;
        }

        //found a prime!  iterate through all products that start with this one
        //They can only contain the same or smaller primes

        //current product
        uint32_t curprod = (uint32_t)test;
        indexes[0] = nprimes;
        products[0] = curprod;
        stksize = 1;

        //remember the found prime (first time through, nprimes == NUMPREMARK)
        primes[nprimes++] = curprod;

        //when we extend the product, we add the first non-premarked prime
        uint32_t extensionPrime = primes[NUMPREMARK];
        //which we can only do if the current product is at most this big
        uint32_t extensionMax = MAXTEST/primes[NUMPREMARK];

        while (curprod <= extensionMax)
        {
            //extend the product with the smallest non-premarked prime
            indexes[stksize]=NUMPREMARK;
            products[stksize++]=(curprod*=extensionPrime);
            SETBIT(mem, curprod);
        }

        for (;;)
        {
            //Can't extend current product.
            //Pop the stacks until we get to a factor we can increase while keeping
            //the factors in descending order and keeping the product small enough
            if (--stksize <= 0)
            {
                //didn't find one
                break;
            }
            if (indexes[stksize]>=indexes[stksize-1])
            {
                //can't increase this factor without breaking descending order
                continue;
            }

            uint64_t testprod=products[stksize-1];
            testprod*=primes[++(indexes[stksize])];
            if (testprod>MAXTEST)
            {
                //can't increase this factor without exceeding our array size
                continue;
            }
            //yay! - we can increment here to the next composite product
            curprod=(uint32_t)testprod;
            products[stksize++] = curprod;
            SETBIT(mem, curprod);

            while (curprod <= extensionMax)
            {
                //extend the product with the smallest non-premarked prime
                indexes[stksize]=NUMPREMARK;
                products[stksize++]=(curprod*=extensionPrime);
                SETBIT(mem, curprod);
            }
        }
    }
    clock_t end_t = clock();
    printf("Found %ld primes\n", nprimes);
    free(mem);
    free(primes);
    printf("Time: %f\n", (double)(end_t - start_t) / CLOCKS_PER_SEC);
}

请注意，我的筛子是从筛子或埃拉托色尼筛子开始的，比你的筛子优化得好一点。主要区别在于我们只为位掩码数组中的奇数分配位。该部分的速度差异并不显着。