c++ - 将优化的埃拉托色尼筛从 Python 移植到 C++

Question

不久前，我在 python 中使用了(速度极快的)primesieve，我在这里找到了它:Fastest way to list all primes below 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 的倍数”，然后用 1000 行 C 文件猛击你)会很棒。

我认为 100% 直接和文字移植不会有问题，但毕竟这是为了学习，这完全没用。

编辑

看过原始 numpy 版本中的代码后，它实际上很容易实现，并且经过一些思考也不太难理解。这是我想出的C++版本。我将其完整版本发布在这里，以帮助更多的读者，以防他们需要一个不超过 200 万行代码的非常高效的 p​​rimesieve。这个 primesieve 在与上述相同的机器上大约 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);
    }
}

Answer 1

我会尽力解释。 sieve 数组有一个不寻常的索引方案；它为每个与 1 或 5 mod 6 全等的数字存储一个位。因此，数字 6*k + 1 将存储在位置 2*k 和 k*6 + 5 将存储在位置 2*k + 1 中。 3*i+1|1 运算与此相反:它采用 2*n 形式的数字并将其转换为 6*n + 1 ，并采用 2*n + 1 并将其转换为 6*n + 5 (+1|1 事物转换0 到 1 和 3 到 5)。主循环迭代 k 遍历具有该属性的所有数字，从 5 开始(当 i 为 1 时)； i 是 sieve 中数字 k 的相应索引。第一个切片更新为 sieve，然后清除筛子中索引形式为 k*k/3 + 2*m*k 的所有位(对于 m 自然数)；这些索引的相应数字从 k^2 开始，并在每一步增加 6*k。第二个切片更新从索引 k*(k-2*(i&1)+4)/3 开始(对于 k，数字 k * (k+4) 与 1 mod 6 一致，否则 k * (k+2))，同样将数字增加 6 *k 每一步。

这是另一种解释尝试:让candidates是至少为5并且与1或5全等的所有数字的集合 mod 6。如果将该集合中的两个元素相乘，则会得到该集合中的另一个元素。令 candidates 中某些 k 的 succ(k) 为 candidates 中的下一个元素(按数字顺序)大于k。在这种情况下，筛子的内部循环基本上是(使用 sieve 的正常索引):

for k in candidates:
  for (l = k; ; l += 6) sieve[k * l] = False
  for (l = succ(k); ; l += 6) sieve[k * l] = False

由于sieve中存储元素的限制，这与:

for k in candidates:
  for l in candidates where l >= k:
    sieve[k * l] = False

它将在某个时刻从筛子中删除候选中k的所有倍数(除了k本身)(无论当前k 之前用作 l 或现在用作 k 时)。