将优化的 Eratosthenes 筛从 Python 移植到 C++答案

【问题标题】：Porting optimized Sieve of Eratosthenes from Python to C++将优化的 Eratosthenes 筛从 Python 移植到 C++
【发布时间】：2011-03-13 23:27:01
【问题描述】：

前段时间我在 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 实现的作用，或者您可以为我移植它（虽然帮助不大）。

编辑

一些关于我觉得困难的额外信息。

我对校正变量等所使用的技术以及它们如何组合在一起存在疑问。一个指向解释不同 Eratosthenes 优化的网站的链接（除了简单的网站说“你只需跳过 2、3 和 5 的倍数”，然后用 1000 行 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://*.com/questions/2068372/fastest-way-to-list-all-primes-below-n-in-python/3035188#3035188
// https://*.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);
    }
}

【问题讨论】：

Python 代码没有将单个元素推入数组；它使用一组元素对其进行初始化，然后就地修改它们。这在 C++ 中也会快得多。像你一样从dynamic_bitset 开始，只需填写它，然后在一个循环中写出筛子。即使不尝试从 Python 版本复制技术，这也可能会使您的版本运行得更快。
您在理解 Python 实现的哪个部分时遇到了问题？请将此信息添加到问题中。此外，人们真的不应该为你回答端口，因为就像你说的那样，这对你的理解没有多大帮助。
@Jeremiah Willcock：不正确。在筛分过程中，它正在就地修改，但最终它将单个元素推入数组（列表）return [2,3] + [3*i+1|1 for i in xrange(1,n/3-correction) if sieve[i]]。这几乎就是我正在做的事情。第一个循环是在适当的位置进行筛选，并且已经推送了我们知道的素数，第二个循环是推送在 sievelimit 的根之后的素数。
@Merlyn Morgan-Graham：我对校正变量等所使用的技术以及它们如何组合在一起遇到了麻烦。一个指向解释不同 Eratosthenes 优化的网站的链接（除了简单的“你只需跳过 2、3 和 5 的倍数”然后被 1000 行 C 文件猛烈抨击）会很棒。（编辑：只需阅读“请将此添加到问题中”，现在就这样做。）
@nightcracker：我不认为sieve 的更新本身会改变它的长度，尽管我不能确定。但是，正如您所演示的，它们确实会在返回时添加元素。

标签： c++ python sieve-of-eratosthenes

【解决方案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+4) 的 k 与 1 一致，否则为 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

这将在某个时刻从筛子中删除candidates 中k 的所有倍数（k 本身除外）（当当前k 更早用作l 时或使用它时）现在就像k）。

【讨论】：

好吧，这基本上就是我为我的筛子所做的排除 2 的操作，但接下来会稍微复杂一些。我会通读几遍，看看是否可以将其转换为 C++ 实现。在每种情况下都为努力 +1。
@nightcracker：我现在完成了索引->数字计算。代码很棘手，但这主要是因为sieve 上的索引。
@nightcracker：看看我的第二部分，如果这更有意义。

【解决方案2】：

紧随 Howard Hinnant 的回答，霍华德，您不必测试不能被 2、3 或 5 整除的所有自然数集合中的数字本身是否为素数。您只需将数组中的每个数字（自消除的 1 除外）乘以自身和数组中的每个后续数字。这些重叠的产品将为您提供数组中的所有非素数，直到您扩展确定性乘法过程的任何点。因此，数组中的第一个非素数将是 7 平方或 49。第二个、7 乘以 11 或 77 等。这里有完整的解释：http://www.primesdemystified.com

【讨论】：

【解决方案3】：

顺便说一句，您可以“近似”质数。调用近似素数P。这里有几个公式：

P = 2*k+1 // 不能被 2 整除

P = 6*k + {1, 5} // 不能整除 2, 3

P = 30*k + {1, 7, 11, 13, 17, 19, 23, 29} // 不能被 2, 3, 5 整除

通过这些公式找到的数字集合的性质是 P 可能不是素数，但是所有素数都在集合 P 中。如果你只测试集合 P 中的素数，你不会错过任何一个。

您可以将这些公式重新表述为：

P = X*k + {-i, -j, -k, k, j, i}

如果这样对你更方便的话。

Here 是一些使用这种技术的代码，其中 P 的公式不能被 2、3、5、7 整除。

此链接可能代表该技术的实际应用程度。

【讨论】：