首页
学习
活动
专区
工具
TVP
发布
社区首页 >问答首页 >用递归公式改进纯Python素数筛

用递归公式改进纯Python素数筛
EN

Stack Overflow用户
提问于 2010-07-20 06:01:27
回答 2查看 3.1K关注 0票数 16

我试图通过去掉子列表长度的复杂公式来进一步优化素数线程中的冠军解决方案。相同子序列的len()太慢,因为len的开销很大,而且生成子序列的开销也很大。这看起来会稍微提高函数的速度,但我还不能去掉除法,尽管我只在条件语句中做了除法。当然,我可以尝试通过对n而不是n*n进行起始标记的优化来简化长度计算。

我将除法/替换为整数除法//以兼容Python 3或

代码语言:javascript
复制
from __future__ import division

如果这个递归公式可以帮助加速numpy解,我也会很感兴趣,但我没有太多使用numpy的经验。

如果您为代码启用了special,那么故事就会完全不同,而且Atkins sieve代码会比这种特殊的切片技术更快。

代码语言:javascript
复制
import cProfile

def rwh_primes1(n):
    # http://stackoverflow.com/questions/2068372/fastest-way-to-list-all-primes-below-n-in-python/3035188#3035188
    """ Returns  a list of primes < n """
    sieve = [True] * (n//2)
    for i in xrange(3,int(n**0.5)+1,2):
        if sieve[i//2]:
            sieve[i*i//2::i] = [False] * ((n-i*i-1)//(2*i)+1)
    return [2] + [2*i+1 for i in xrange(1,n/2) if sieve[i]]

def primes(n):
    # http://stackoverflow.com/questions/2068372/fastest-way-to-list-all-primes-below-n-in-python/3035188#3035188
    # recurrence formula for length by amount1 and amount2 Tony Veijalainen 2010
    """ Returns  a list of primes < n """
    sieve = [True] * (n//2)
    amount1 = n-10
    amount2 = 6

    for i in xrange(3,int(n**0.5)+1,2):
        if sieve[i//2]:
             ## can you make recurrence formula for whole reciprocal?
            sieve[i*i//2::i] = [False] * (amount1//amount2+1)
        amount1-=4*i+4
        amount2+=4

    return [2] + [2*i+1 for i in xrange(1,n//2) if sieve[i]]

numprimes=1000000
print('Profiling')
cProfile.Profile.bias = 4e-6
for test in (rwh_primes1, primes):
    cProfile.run("test(numprimes)")

分析(不同版本之间差别不大)

代码语言:javascript
复制
         3 function calls in 0.191 CPU seconds

   Ordered by: standard name

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
        1    0.006    0.006    0.191    0.191 <string>:1(<module>)
        1    0.185    0.185    0.185    0.185 myprimes.py:3(rwh_primes1)
        1    0.000    0.000    0.000    0.000 {method 'disable' of '_lsprof.Profiler' objects}


         3 function calls in 0.192 CPU seconds

   Ordered by: standard name

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
        1    0.006    0.006    0.192    0.192 <string>:1(<module>)
        1    0.186    0.186    0.186    0.186 myprimes.py:12(primes)
        1    0.000    0.000    0.000    0.000 {method 'disable' of '_lsprof.Profiler' objects}

有趣的是,通过将限制增加到10**8,并将计时修饰器添加到删除分析的函数中:

代码语言:javascript
复制
rwh_primes1 took 23.670 s
primes took 22.792 s
primesieve took 10.850 s

有趣的是,如果你不生成素数列表,而是返回筛子本身,时间大约是数字列表版本的一半。

EN

回答 2

Stack Overflow用户

发布于 2011-04-13 17:07:19

你可以做一个轮子优化。2和3的倍数不是素数,所以根本不要存储它们。然后,您可以从5开始,通过以2、4、2、4、2、4等步长递增来跳过2和3的倍数。

下面是它的C++代码。希望这能有所帮助。

代码语言:javascript
复制
void sieve23()
{
    int lim=sqrt(MAX);
    for(int i=5,bit1=0;i<=lim;i+=(bit1?4:2),bit1^=1)
    {
        if(!isComp[i/3])
        {
            for(int j=i,bit2=1;;)
            {
                j+=(bit2?4*i:2*i);
                bit2=!bit2;
                if(j>=MAX)break;
                isComp[j/3]=1;
            }
        }
    }
}
票数 1
EN

Stack Overflow用户

发布于 2011-04-13 17:33:47

如果您可能决定使用C++来提高速度,我将Python sieve移植到了C++。完整的讨论可以在这里找到:Porting optimized Sieve of Eratosthenes from Python to C++

在使用g++ -O3和N=100000000编译的英特尔Q6600 Ubuntu 10.10上,这需要415毫秒。

代码语言:javascript
复制
#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, sievemax, sievemaxroot;

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

    sievemaxroot = isqrt(N)/3;

    boost::dynamic_bitset<> sieve(sievemax);
    sieve.set();
    sieve[0] = 0;

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

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

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

}
票数 0
EN
页面原文内容由Stack Overflow提供。腾讯云小微IT领域专用引擎提供翻译支持
原文链接:

https://stackoverflow.com/questions/3285443

复制
相关文章

相似问题

领券
问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档