2

前段时间我在 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=100000000g++ -O31230 毫秒。

现在,我希望能在理解上述 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://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);
    }
}
4

3 回答 3

3

我会尽力解释。sieve数组有一个不寻常的索引方案;它为每个等于 1 或 5 mod 6 的数字存储一个位。因此,一个数字6*k + 1将存储在 position2*k并将k*6 + 5存储在 position 2*k + 1。该3*i+1|1操作与此相反:它获取表单的数字2*n并将它们转换为6*n + 1,然后获取2*n + 1并将其转换为(6*n + 5事物+1|1转换01和)。主循环遍历具有该属性的所有数字,从(when is 1);开始 是数字的对应索引。第一个切片更新为35k5iisieveksieve然后用以下形式的索引清除筛子中的所有位k*k/3 + 2*m*k(对于m自然数);这些索引的相应数字从每一步开始k^2并增加6*k。第二个切片更新从索引k*(k-2*(i&1)+4)/3(与 mod 一致的数字等)开始k * (k+4),并且在每个步骤中类似地增加数字。k16k * (k+2)6*k

这是解释的另一种尝试:让candidates是至少为 5 并且与 any15mod一致的所有数字的集合6。如果将该集合中的两个元素相乘,则会得到该集合中的另一个元素。让succ(k)for some kin是大于candidates的下一个元素(按数字顺序)。在这种情况下,筛子的内循环基本上是(使用普通索引 for ):candidatesksieve

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

这将在某个时候从筛子中删除所有倍数的kin candidates(除了k它本身)(无论是当前k使用l之前还是k现在使用时)。

于 2011-03-14T00:20:44.120 回答
1

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

于 2011-03-22T05:40:30.333 回答
0

顺便说一句,您可以“近似”素数。调用近似素数 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中的数字是否为素数,那么所有素数都在集合PIe中,您不会错过任何一个。

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

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

如果这对你更方便。

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

这个链接可能代表了这种技术在多大程度上可以被实际利用。

于 2011-03-14T00:34:31.427 回答