前段时间,我在 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% 左右)。
比较结果(希望这个问题能够成功回答),当前执行时间与N=100000000
, g++ -O3
在 Q6600 Ubuntu 10.10 上为 1230 毫秒。
现在我希望得到一些帮助,帮助我理解上面的 Python 实现的作用,或者你可以帮我移植它(虽然没有那么有帮助)。
EDIT
关于我觉得困难的一些额外信息。
我对校正变量等所使用的技术以及一般情况下如何组合在一起感到困难。一个解释不同埃拉托色尼优化的网站的链接(除了那些简单的网站说“好吧,你只需跳过 2、3 和 5 的倍数”,然后用 1000 行 C 文件猛击你)会很棒。
我不认为我会对 100% 直接和文字移植有任何问题,但毕竟这是为了学习,这是完全没有用的。
EDIT
查看原始 numpy 版本中的代码后,它实际上非常容易实现,并且经过一些思考也不太难理解。这是我想出的C++版本。我将其完整版本发布在这里,以帮助更多的读者,以防他们需要一个不超过 200 万行代码的非常高效的 primesieve。这个 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);
}
}