将优化的埃拉托斯特尼筛法从Python移植到C++

Question

不久前，我在Python中使用了一个非常快速的质数筛选工具，叫做primesieve，具体可以在这里找到：列出所有小于N的质数的最快方法

为了更准确地说，这个实现是：

def primes2(n):
    """ Input n>=6, Returns a list of primes, 2 <= p < n """
    n, correction = n-n%6+6, 2-(n%6>1)
    sieve = [True] * (n/3)
    for i in xrange(1,int(n**0.5)/3+1):
      if sieve[i]:
        k=3*i+1|1
        sieve[      k*k/3      ::2*k] = [False] * ((n/6-k*k/6-1)/k+1)
        sieve[k*(k-2*(i&1)+4)/3::2*k] = [False] * ((n/6-k*(k-2*(i&1)+4)/6-1)/k+1)
    return [2,3] + [3*i+1|1 for i in xrange(1,n/3-correction) if sieve[i]]

现在我大致明白了通过自动跳过2、3等的倍数来优化的思路，但当我尝试把这个算法移植到C++时就遇到了困难（我对Python理解得不错，对C++的理解一般/较差，但也足够用）。

我目前自己写的代码是这样的（isqrt()只是一个简单的整数平方根函数）：

template <class T>
void primesbelow(T N, std::vector<T> &primes) {
    T sievemax = (N-3 + (1-(N % 2))) / 2;
    T i;
    T sievemaxroot = isqrt(sievemax) + 1;

    boost::dynamic_bitset<> sieve(sievemax);
    sieve.set();

    primes.push_back(2);

    for (i = 0; i <= sievemaxroot; i++) {
        if (sieve[i]) {
            primes.push_back(2*i+3);
            for (T j = 3*i+3; j <= sievemax; j += 2*i+3) sieve[j] = 0; // filter multiples
        }
    }

    for (; i <= sievemax; i++) {
        if (sieve[i]) primes.push_back(2*i+3);
    }
}

这个实现还不错，自动跳过了2的倍数，但如果我能把Python的实现移植过来，我觉得速度会快很多（大约快50%-30%）。

为了比较结果（希望这个问题能得到成功的解答），在一台Q6600的Ubuntu 10.10上，使用N=100000000和g++ -O3的当前执行时间是1230毫秒。

我现在希望能得到一些帮助，要么理解上面Python实现的具体做法，要么帮我移植一下（不过后者的帮助性不大）。

编辑

关于我觉得困难的地方，补充一些信息。

我对使用的技术，比如修正变量，以及整体是如何结合在一起的，感到困惑。如果有链接能解释不同的埃拉托斯特尼筛法优化（除了那些简单说“你只需跳过2、3和5的倍数”然后给你一千行C代码的网站）就太好了。

我觉得如果是100%直接和字面上的移植应该不会有问题，但毕竟这是为了学习，那样就完全没用。

编辑

在查看了原始numpy版本的代码后，实际上实现起来相当简单，经过一些思考也不难理解。这是我想到的C++版本。我把它完整地发布在这里，以帮助后来的读者，万一他们需要一个效率不错的质数筛选工具，而不是两百万行的代码。这个质数筛选工具在上面提到的同一台机器上，能在大约415毫秒内找出所有小于100000000的质数。这是3倍的速度提升，比我预期的要好！

#include <vector>
#include <boost/dynamic_bitset.hpp>

// http://vault.embedded.com/98/9802fe2.htm - integer square root
unsigned short isqrt(unsigned long a) {
    unsigned long rem = 0;
    unsigned long root = 0;

    for (short i = 0; i < 16; i++) {
        root <<= 1;
        rem = ((rem << 2) + (a >> 30));
        a <<= 2;
        root++;

        if (root <= rem) {
            rem -= root;
            root++;
        } else root--;

    }

    return static_cast<unsigned short> (root >> 1);
}

// https://stackoverflow.com/questions/2068372/fastest-way-to-list-all-primes-below-n-in-python/3035188#3035188
// https://stackoverflow.com/questions/5293238/porting-optimized-sieve-of-eratosthenes-from-python-to-c/5293492
template <class T>
void primesbelow(T N, std::vector<T> &primes) {
    T i, j, k, l, sievemax, sievemaxroot;

    sievemax = N/3;
    if ((N % 6) == 2) sievemax++;

    sievemaxroot = isqrt(N)/3;

    boost::dynamic_bitset<> sieve(sievemax);
    sieve.set();

    primes.push_back(2);
    primes.push_back(3);

    for (i = 1; i <= sievemaxroot; i++) {
        if (sieve[i]) {
            k = (3*i + 1) | 1;
            l = (4*k-2*k*(i&1)) / 3;

            for (j = k*k/3; j < sievemax; j += 2*k) {
                sieve[j] = 0;
                sieve[j+l] = 0;
            }

            primes.push_back(k);
        }
    }

    for (i = sievemaxroot + 1; i < sievemax; i++) {
        if (sieve[i]) primes.push_back((3*i+1)|1);
    }
}

性能提升编程学习代码效率质数筛选埃拉托斯特尼筛法算法移植 C++优化整数平方根

将优化的埃拉托斯特尼筛法从Python移植到C++

3 个回答

撰写回答