python - 通过递归公式改进纯Python素筛

Question

我正在尝试通过取出子列表长度的复杂公式来进一步优化素数线程中的冠军解决方案。同一子序列的 len() 太慢，因为 len 很昂贵并且生成子序列很昂贵。这看起来稍微加快了函数的速度，但我还不能取消除法，尽管我只在条件语句中进行除法。当然，我可以尝试通过取消开始标记为 n 而不是 n*n 的优化来简化长度计算......

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

兼容

from __future__ import division

另外，如果这个递归公式可以帮助加快 numpy 的解决方案，我会很感兴趣，但我没有太多使用 numpy 的经验。

如果为代码启用 psyco，情况将完全不同，但是阿特金斯筛代码变得比这种特殊的切片技术更快。

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)")

Profiling(版本之间差别不大)

         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 并将计时装饰器添加到函数中以移除分析:

rwh_primes1 took 23.670 s
primes took 22.792 s
primesieve took 10.850 s

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

Answer 1

您可以进行车轮优化。 2 和 3 的倍数不是质数，所以根本不要存储它们。然后，您可以从 5 开始，通过以 2、4、2、4、2、4 等步长递增来跳过 2 和 3 的倍数。

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

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;
            }
        }
    }
}