40

我使用埃拉托色尼筛法和 Python 3.1编写了一个素数生成器。该代码在ideone.com上以 0.32 秒的时间正确而优雅地运行,以生成高达 1,000,000 的素数。

# from bitstring import BitString

def prime_numbers(limit=1000000):
    '''Prime number generator. Yields the series
    2, 3, 5, 7, 11, 13, 17, 19, 23, 29 ...    
    using Sieve of Eratosthenes.
    '''
    yield 2
    sub_limit = int(limit**0.5) 
    flags = [False, False] + [True] * (limit - 2)   
#     flags = BitString(limit)
    # Step through all the odd numbers
    for i in range(3, limit, 2):       
        if flags[i] is False:
#        if flags[i] is True:
            continue
        yield i
        # Exclude further multiples of the current prime number
        if i <= sub_limit:
            for j in range(i*3, limit, i<<1):
                flags[j] = False
#                flags[j] = True

问题是,当我尝试生成高达 1,000,000,000 的数字时,内存不足。

    flags = [False, False] + [True] * (limit - 2)   
MemoryError

可以想象,分配 10 亿个布尔值( Python 中每个1 字节4 或 8 字节(见注释))确实不可行,所以我研究了bitstring。我想,为每个标志使用 1 位会更节省内存。然而,该程序的性能急剧下降 - 运行时间为 24 秒,质数高达 1,000,000。这可能是由于位串的内部实现。

您可以注释/取消注释这三行,以查看我更改为使用 BitString 的内容,如上面的代码片段。

我的问题是,有没有办法加快我的程序,有或没有位串?

编辑:请在发布之前自己测试代码。自然,我不能接受运行速度比我现有代码慢的答案。

再次编辑:

我已经在我的机器上编译了一个基准测试列表。

4

11 回答 11

34

您的版本有几个小的优化。通过颠倒 True 和 False 的角色,您可以将“ if flags[i] is False:”更改为“ if flags[i]:”。并且第二条语句的起始值range可以i*i代替i*3. 您的原始版本在我的系统上需要 0.166 秒。通过这些更改,以下版本在我的系统上需要 0.156 秒。

def prime_numbers(limit=1000000):
    '''Prime number generator. Yields the series
    2, 3, 5, 7, 11, 13, 17, 19, 23, 29 ...
    using Sieve of Eratosthenes.
    '''
    yield 2
    sub_limit = int(limit**0.5)
    flags = [True, True] + [False] * (limit - 2)
    # Step through all the odd numbers
    for i in range(3, limit, 2):
        if flags[i]:
            continue
        yield i
        # Exclude further multiples of the current prime number
        if i <= sub_limit:
            for j in range(i*i, limit, i<<1):
                flags[j] = True

不过,这对您的记忆问题没有帮助。

进入 C 扩展的世界,我使用了gmpy的开发版本。(免责声明:我是维护者之一。)开发版本称为 gmpy2 并支持称为 xmpz 的可变整数。使用 gmpy2 和以下代码,我的运行时间为 0.140 秒。限制为 1,000,000,000 的运行时间为 158 秒。

import gmpy2

def prime_numbers(limit=1000000):
    '''Prime number generator. Yields the series
    2, 3, 5, 7, 11, 13, 17, 19, 23, 29 ...
    using Sieve of Eratosthenes.
    '''
    yield 2
    sub_limit = int(limit**0.5)
    # Actual number is 2*bit_position + 1.
    oddnums = gmpy2.xmpz(1)
    current = 0
    while True:
        current += 1
        current = oddnums.bit_scan0(current)
        prime = 2 * current + 1
        if prime > limit:
            break
        yield prime
        # Exclude further multiples of the current prime number
        if prime <= sub_limit:
            for j in range(2*current*(current+1), limit>>1, prime):
                oddnums.bit_set(j)

推动优化并牺牲清晰度,我使用以下代码获得了 0.107 和 123 秒的运行时间:

import gmpy2

def prime_numbers(limit=1000000):
    '''Prime number generator. Yields the series
    2, 3, 5, 7, 11, 13, 17, 19, 23, 29 ...
    using Sieve of Eratosthenes.
    '''
    yield 2
    sub_limit = int(limit**0.5)
    # Actual number is 2*bit_position + 1.
    oddnums = gmpy2.xmpz(1)
    f_set = oddnums.bit_set
    f_scan0 = oddnums.bit_scan0
    current = 0
    while True:
        current += 1
        current = f_scan0(current)
        prime = 2 * current + 1
        if prime > limit:
            break
        yield prime
        # Exclude further multiples of the current prime number
        if prime <= sub_limit:
            list(map(f_set,range(2*current*(current+1), limit>>1, prime)))

编辑:基于这个练习,我将 gmpy2 修改为 accept xmpz.bit_set(iterator)。使用以下代码,所有小于 1,000,000,000 的素数的运行时间对于 Python 2.7 为 56 秒,对于 Python 3.2 为 74 秒。(如评论中所述,xrange比 快range。)

import gmpy2

try:
    range = xrange
except NameError:
    pass

def prime_numbers(limit=1000000):
    '''Prime number generator. Yields the series
    2, 3, 5, 7, 11, 13, 17, 19, 23, 29 ...
    using Sieve of Eratosthenes.
    '''
    yield 2
    sub_limit = int(limit**0.5)
    oddnums = gmpy2.xmpz(1)
    f_scan0 = oddnums.bit_scan0
    current = 0
    while True:
        current += 1
        current = f_scan0(current)
        prime = 2 * current + 1
        if prime > limit:
            break
        yield prime
        if prime <= sub_limit:
            oddnums.bit_set(iter(range(2*current*(current+1), limit>>1, prime)))

编辑#2:再试一次!我修改了 gmpy2 以接受xmpz.bit_set(slice). 使用以下代码,对于 Python 2.7 和 Python 3.2,所有小于 1,000,000,000 的素数的运行时间约为 40 秒。

from __future__ import print_function
import time
import gmpy2

def prime_numbers(limit=1000000):
    '''Prime number generator. Yields the series
    2, 3, 5, 7, 11, 13, 17, 19, 23, 29 ...
    using Sieve of Eratosthenes.
    '''
    yield 2
    sub_limit = int(limit**0.5)
    flags = gmpy2.xmpz(1)
    # pre-allocate the total length
    flags.bit_set((limit>>1)+1)
    f_scan0 = flags.bit_scan0
    current = 0
    while True:
        current += 1
        current = f_scan0(current)
        prime = 2 * current + 1
        if prime > limit:
            break
        yield prime
        if prime <= sub_limit:
            flags.bit_set(slice(2*current*(current+1), limit>>1, prime))

start = time.time()
result = list(prime_numbers(1000000000))
print(time.time() - start)

编辑#3:我已更新 gmpy2 以正确支持 xmpz 位级别的切片。性能没有变化,但 API 非常好。我做了一些调整,我把时间缩短到了大约 37 秒。(请参阅编辑 #4 以了解 gmpy2 2.0.0b1 中的更改。)

from __future__ import print_function
import time
import gmpy2

def prime_numbers(limit=1000000):
    '''Prime number generator. Yields the series
    2, 3, 5, 7, 11, 13, 17, 19, 23, 29 ...
    using Sieve of Eratosthenes.
    '''
    sub_limit = int(limit**0.5)
    flags = gmpy2.xmpz(1)
    flags[(limit>>1)+1] = True
    f_scan0 = flags.bit_scan0
    current = 0
    prime = 2
    while prime <= sub_limit:
        yield prime
        current += 1
        current = f_scan0(current)
        prime = 2 * current + 1
        flags[2*current*(current+1):limit>>1:prime] = True
    while prime <= limit:
        yield prime
        current += 1
        current = f_scan0(current)
        prime = 2 * current + 1

start = time.time()
result = list(prime_numbers(1000000000))
print(time.time() - start)

编辑#4:我在 gmpy2 2.0.0b1 中做了一些改变,打破了前面的例子。gmpy2 不再将 True 视为提供无限源的 1 位的特殊值。-1 应改为使用。

from __future__ import print_function
import time
import gmpy2

def prime_numbers(limit=1000000):
    '''Prime number generator. Yields the series
    2, 3, 5, 7, 11, 13, 17, 19, 23, 29 ...
    using Sieve of Eratosthenes.
    '''
    sub_limit = int(limit**0.5)
    flags = gmpy2.xmpz(1)
    flags[(limit>>1)+1] = 1
    f_scan0 = flags.bit_scan0
    current = 0
    prime = 2
    while prime <= sub_limit:
        yield prime
        current += 1
        current = f_scan0(current)
        prime = 2 * current + 1
        flags[2*current*(current+1):limit>>1:prime] = -1
    while prime <= limit:
        yield prime
        current += 1
        current = f_scan0(current)
        prime = 2 * current + 1

start = time.time()
result = list(prime_numbers(1000000000))
print(time.time() - start)

编辑#5:我对 gmpy2 2.0.0b2 做了一些改进。您现在可以遍历所有设置或清除的位。运行时间提高了约 30%。

from __future__ import print_function
import time
import gmpy2

def sieve(limit=1000000):
    '''Returns a generator that yields the prime numbers up to limit.'''

    # Increment by 1 to account for the fact that slices do not include
    # the last index value but we do want to include the last value for
    # calculating a list of primes.
    sieve_limit = gmpy2.isqrt(limit) + 1
    limit += 1

    # Mark bit positions 0 and 1 as not prime.
    bitmap = gmpy2.xmpz(3)

    # Process 2 separately. This allows us to use p+p for the step size
    # when sieving the remaining primes.
    bitmap[4 : limit : 2] = -1

    # Sieve the remaining primes.
    for p in bitmap.iter_clear(3, sieve_limit):
        bitmap[p*p : limit : p+p] = -1

    return bitmap.iter_clear(2, limit)

if __name__ == "__main__":
    start = time.time()
    result = list(sieve(1000000000))
    print(time.time() - start)
    print(len(result))
于 2010-05-25T03:41:49.593 回答
10

好的,这是我的第二个答案,但由于速度至关重要,我认为我不得不提到bitarray模块——即使它是bitstring的克星:)。它非常适合这种情况,因为它不仅是一个 C 扩展(并且比纯 Python 希望更快),而且它还支持切片分配。不过,它还不适用于 Python 3。

我什至没有尝试优化这个,我只是重写了位串版本。在我的机器上,对于一百万以下的素数,我得到 0.16 秒。

对于十亿,它运行得非常好,并在 2 分 31 秒内完成。

import bitarray

def prime_bitarray(limit=1000000):
    yield 2
    flags = bitarray.bitarray(limit)
    flags.setall(False)
    sub_limit = int(limit**0.5)
    for i in xrange(3, limit, 2):
        if not flags[i]:
            yield i
            if i <= sub_limit:
                flags[3*i:limit:i*2] = True
于 2010-05-24T20:53:33.490 回答
8

好的,这是我今晚完成的(接近完整的)综合基准测试,以查看哪些代码运行速度最快。希望有人会发现这个列表很有用。我省略了在我的机器上需要超过 30 秒才能完成的任何事情。

我要感谢所有提出意见的人。我从你的努力中获得了很多见识,我希望你也一样。

我的机器:AMD ZM-86,2.40 Ghz 双核,4GB RAM。这是一台 HP Touchsmart Tx2 笔记本电脑。请注意,虽然我可能已链接到 pastebin,但我在自己的机器上对以下所有内容进行了基准测试。

一旦我能够构建它,我将添加 gmpy2 基准。

所有基准测试均在 Python 2.6 x86 中测试

返回一个质数列表 n 最多 1,000,000:(使用Python 生成器)

Sebastian 的 numpy 生成器版本(更新) - 121 ms @

马克筛 + 轮- 154 毫秒

带切片的罗伯特版本- 159 毫秒

我的改进版切片 - 205 毫秒

带枚举的 Numpy 生成器- 249 ms @

马克的基本筛- 317 毫秒

casevh 对我原来的解决方案的改进- 343 毫秒

我修改后的 numpy 生成器解决方案- 407 ms

我在问题中的原始方法- 409 ms

Bitarray 解决方案- 414 毫秒 @

带字节数组的纯 Python - 1394 毫秒 @

Scott 的 BitString 解决方案- 6659 毫秒 @

'@' 表示此方法能够在我的机器设置上生成最多 n < 1,000,000,000。

此外,如果您不需要生成器并且只想要一次整个列表:

RosettaCode 的 numpy 解决方案- 32 ms @

(numpy 解决方案还能够生成多达 10 亿个,这需要 61.6259 秒。我怀疑内存被交换了一次,因此是双倍的时间。)

于 2010-05-25T10:22:40.187 回答
6

相关问题:在 python 中列出 N 以下所有素数的最快方法

嗨,我也在寻找 Python 中的代码以尽可能快地生成高达10^9的素数,由于内存问题,这很困难。到目前为止,我想出这个来生成高达10^610^7的素数(在我的旧机器上分别计时 0,171s 和 1,764s),这似乎比“我的上一篇文章中带有切片的改进版本,可能是因为我使用n//i-i +1(见下文)而不是其他代码中的类似len(flags[i2::i<<1])代码。如果我错了,请纠正我。非常欢迎任何改进建议。

def primes(n):
    """ Returns  a list of primes < n """
    sieve = [True] * n
    for i in xrange(3,int(n**0.5)+1,2):
        if sieve[i]:
            sieve[i*i::2*i]=[False]*((n-i*i-1)/(2*i)+1)
    return [2] + [i for i in xrange(3,n,2) if sieve[i]]

在他的一个代码中,Xavier 使用flags[i2::i<<1]len(flags[i2::i<<1])i*i..n我明确地计算了大小,使用我们之间有(n-i*i)//2*i倍数的事实,2*i因为我们也想计算i*i我们总和1所以len(sieve[i*i::2*i])等于(n-i*i)//(2*i) +1。这使代码更快。基于上述代码的基本生成器将是:

def primesgen(n):
    """ Generates all primes <= n """
    sieve = [True] * n
    yield 2
    for i in xrange(3,int(n**0.5)+1,2):
        if sieve[i]:
            yield i
            sieve[i*i::2*i] = [False]*((n-i*i-1)/(2*i)+1)
    for i in xrange(i+2,n,2):
        if sieve[i]: yield i

稍加修改,就可以编写上面代码的稍慢版本,从一半大小的筛子开始,sieve = [True] * (n//2)并且工作相同n。不确定这将如何反映在内存问题中。作为实现示例,这里是 numpy Rosetta 代码的修改版本(可能更快),从一半大小的筛子开始。

import numpy
def primesfrom3to(n):
    """ Returns a array of primes, 3 <= p < n """
    sieve = numpy.ones(n/2, dtype=numpy.bool)
    for i in xrange(3,int(n**0.5)+1,2):
        if sieve[i/2]: sieve[i*i/2::i] = False
    return 2*numpy.nonzero(sieve)[0][1::]+1

更快、更明智的生成器将是:

import numpy
def primesgen1(n):
""" Input n>=6, Generates all primes < n """
sieve1 = numpy.ones(n/6+1, dtype=numpy.bool)
sieve5 = numpy.ones(n/6  , dtype=numpy.bool)
sieve1[0] = False
yield 2; yield 3;
for i in xrange(int(n**0.5)/6+1):
    if sieve1[i]:
        k=6*i+1; yield k;
        sieve1[ ((k*k)/6) ::k] = False
        sieve5[(k*k+4*k)/6::k] = False
    if sieve5[i]:
        k=6*i+5; yield k;
        sieve1[ ((k*k)/6) ::k] = False
        sieve5[(k*k+2*k)/6::k] = False
for i in xrange(i+1,n/6):
        if sieve1[i]: yield 6*i+1
        if sieve5[i]: yield 6*i+5
if n%6 > 1:
    if sieve1[i+1]: yield 6*(i+1)+1

或更多代码:

import numpy
def primesgen(n):
    """ Input n>=30, Generates all primes < n """
    size = n/30 + 1
    sieve01 = numpy.ones(size, dtype=numpy.bool)
    sieve07 = numpy.ones(size, dtype=numpy.bool)
    sieve11 = numpy.ones(size, dtype=numpy.bool)
    sieve13 = numpy.ones(size, dtype=numpy.bool)
    sieve17 = numpy.ones(size, dtype=numpy.bool)
    sieve19 = numpy.ones(size, dtype=numpy.bool)
    sieve23 = numpy.ones(size, dtype=numpy.bool)
    sieve29 = numpy.ones(size, dtype=numpy.bool)
    sieve01[0] = False
    yield 2; yield 3; yield 5;
    for i in xrange(int(n**0.5)/30+1):
        if sieve01[i]:
            k=30*i+1; yield k;
            sieve01[     (k*k)/30::k] = False
            sieve07[(k*k+ 6*k)/30::k] = False
            sieve11[(k*k+10*k)/30::k] = False
            sieve13[(k*k+12*k)/30::k] = False
            sieve17[(k*k+16*k)/30::k] = False
            sieve19[(k*k+18*k)/30::k] = False
            sieve23[(k*k+22*k)/30::k] = False
            sieve29[(k*k+28*k)/30::k] = False
        if sieve07[i]:
            k=30*i+7; yield k;
            sieve01[(k*k+ 6*k)/30::k] = False
            sieve07[(k*k+24*k)/30::k] = False
            sieve11[(k*k+16*k)/30::k] = False
            sieve13[(k*k+12*k)/30::k] = False
            sieve17[(k*k+ 4*k)/30::k] = False
            sieve19[     (k*k)/30::k] = False
            sieve23[(k*k+22*k)/30::k] = False
            sieve29[(k*k+10*k)/30::k] = False
        if sieve11[i]:
            k=30*i+11; yield k;
            sieve01[     (k*k)/30::k] = False
            sieve07[(k*k+ 6*k)/30::k] = False
            sieve11[(k*k+20*k)/30::k] = False
            sieve13[(k*k+12*k)/30::k] = False
            sieve17[(k*k+26*k)/30::k] = False
            sieve19[(k*k+18*k)/30::k] = False
            sieve23[(k*k+ 2*k)/30::k] = False
            sieve29[(k*k+ 8*k)/30::k] = False
        if sieve13[i]:
            k=30*i+13; yield k;
            sieve01[(k*k+24*k)/30::k] = False
            sieve07[(k*k+ 6*k)/30::k] = False
            sieve11[(k*k+ 4*k)/30::k] = False
            sieve13[(k*k+18*k)/30::k] = False
            sieve17[(k*k+16*k)/30::k] = False
            sieve19[     (k*k)/30::k] = False
            sieve23[(k*k+28*k)/30::k] = False
            sieve29[(k*k+10*k)/30::k] = False
        if sieve17[i]:
            k=30*i+17; yield k;
            sieve01[(k*k+ 6*k)/30::k] = False
            sieve07[(k*k+24*k)/30::k] = False
            sieve11[(k*k+26*k)/30::k] = False
            sieve13[(k*k+12*k)/30::k] = False
            sieve17[(k*k+14*k)/30::k] = False
            sieve19[     (k*k)/30::k] = False
            sieve23[(k*k+ 2*k)/30::k] = False
            sieve29[(k*k+20*k)/30::k] = False
        if sieve19[i]:
            k=30*i+19; yield k;
            sieve01[     (k*k)/30::k] = False
            sieve07[(k*k+24*k)/30::k] = False
            sieve11[(k*k+10*k)/30::k] = False
            sieve13[(k*k+18*k)/30::k] = False
            sieve17[(k*k+ 4*k)/30::k] = False
            sieve19[(k*k+12*k)/30::k] = False
            sieve23[(k*k+28*k)/30::k] = False
            sieve29[(k*k+22*k)/30::k] = False
        if sieve23[i]:
            k=30*i+23; yield k;
            sieve01[(k*k+24*k)/30::k] = False
            sieve07[(k*k+ 6*k)/30::k] = False
            sieve11[(k*k+14*k)/30::k] = False
            sieve13[(k*k+18*k)/30::k] = False
            sieve17[(k*k+26*k)/30::k] = False
            sieve19[     (k*k)/30::k] = False
            sieve23[(k*k+ 8*k)/30::k] = False
            sieve29[(k*k+20*k)/30::k] = False
        if sieve29[i]:
            k=30*i+29; yield k;
            sieve01[     (k*k)/30::k] = False
            sieve07[(k*k+24*k)/30::k] = False
            sieve11[(k*k+20*k)/30::k] = False
            sieve13[(k*k+18*k)/30::k] = False
            sieve17[(k*k+14*k)/30::k] = False
            sieve19[(k*k+12*k)/30::k] = False
            sieve23[(k*k+ 8*k)/30::k] = False
            sieve29[(k*k+ 2*k)/30::k] = False
    for i in xrange(i+1,n/30):
            if sieve01[i]: yield 30*i+1
            if sieve07[i]: yield 30*i+7
            if sieve11[i]: yield 30*i+11
            if sieve13[i]: yield 30*i+13
            if sieve17[i]: yield 30*i+17
            if sieve19[i]: yield 30*i+19
            if sieve23[i]: yield 30*i+23
            if sieve29[i]: yield 30*i+29
    i += 1
    mod30 = n%30
    if mod30 > 1:
        if sieve01[i]: yield 30*i+1
    if mod30 > 7:
        if sieve07[i]: yield 30*i+7
    if mod30 > 11:
        if sieve11[i]: yield 30*i+11
    if mod30 > 13:
        if sieve13[i]: yield 30*i+13
    if mod30 > 17:
        if sieve17[i]: yield 30*i+17
    if mod30 > 19:
        if sieve19[i]: yield 30*i+19
    if mod30 > 23:
        if sieve23[i]: yield 30*i+23

Ps:如果您对如何加速此代码有任何建议,从更改操作顺序到预先计算任何内容,请发表评论。

于 2010-05-25T21:02:28.383 回答
4

这是我不久前写的一个版本;与您的速度进行比较可能会很有趣。但是,它对空间问题没有任何作用(实际上,它们可能比您的版本更糟)。

from math import sqrt

def basicSieve(n):
    """Given a positive integer n, generate the primes < n."""
    s = [1]*n
    for p in xrange(2, 1+int(sqrt(n-1))):
        if s[p]:
            a = p*p
            s[a::p] = [0] * -((a-n)//p)
    for p in xrange(2, n):
        if s[p]:
            yield p 

我有更快的版本,使用轮子,但它们要复杂得多。原始来源在这里

好的,这是使用轮子的版本。 wheelSieve是主要的切入点。

from math import sqrt
from bisect import bisect_left

def basicSieve(n):
    """Given a positive integer n, generate the primes < n."""
    s = [1]*n
    for p in xrange(2, 1+int(sqrt(n-1))):
        if s[p]:
            a = p*p
            s[a::p] = [0] * -((a-n)//p)
    for p in xrange(2, n):
        if s[p]:
            yield p

class Wheel(object):
    """Class representing a wheel.

    Attributes:
       primelimit -> wheel covers primes < primelimit.
       For example, given a primelimit of 6
       the wheel primes are 2, 3, and 5.
       primes -> list of primes less than primelimit
       size -> product of the primes in primes;  the modulus of the wheel
       units -> list of units modulo size
       phi -> number of units

    """
    def __init__(self, primelimit):
        self.primelimit = primelimit
        self.primes = list(basicSieve(primelimit))

        # compute the size of the wheel
        size = 1
        for p in self.primes:
            size *= p
        self.size = size

        # find the units by sieving
        units = [1] * self.size
        for p in self.primes:
            units[::p] = [0]*(self.size//p)
        self.units = [i for i in xrange(self.size) if units[i]]

        # number of units
        self.phi = len(self.units)

    def to_index(self, n):
        """Compute alpha(n), where alpha is an order preserving map
        from the set of units modulo self.size to the nonnegative integers.

        If n is not a unit, the index of the first unit greater than n
        is given."""
        return bisect_left(self.units, n % self.size) + n // self.size * self.phi

    def from_index(self, i):
        """Inverse of to_index."""

        return self.units[i % self.phi] + i // self.phi * self.size

def wheelSieveInner(n, wheel):
    """Given a positive integer n and a wheel, find the wheel indices of
    all primes that are less than n, and that are units modulo the
    wheel modulus.
    """

    # renaming to avoid the overhead of attribute lookups
    U = wheel.units
    wS = wheel.size
    # inverse of unit map
    UI = dict((u, i) for i, u in enumerate(U))
    nU = len(wheel.units)

    sqroot = 1+int(sqrt(n-1)) # ceiling of square root of n

    # corresponding index (index of next unit up)
    sqrti = bisect_left(U, sqroot % wS) + sqroot//wS*nU
    upper = bisect_left(U, n % wS) + n//wS*nU
    ind2 = bisect_left(U, 2 % wS) + 2//wS*nU

    s = [1]*upper
    for i in xrange(ind2, sqrti):
        if s[i]:
            q = i//nU
            u = U[i%nU]
            p = q*wS+u
            u2 = u*u
            aq, au = (p+u)*q+u2//wS, u2%wS
            wp = p * nU
            for v in U:
                # eliminate entries corresponding to integers
                # congruent to p*v modulo p*wS
                uvr = u*v%wS
                m = aq + (au > uvr)
                bot = (m + (q*v + u*v//wS - m) % p) * nU + UI[uvr]
                s[bot::wp] = [0]*-((bot-upper)//wp)
    return s

def wheelSieve(n, wheel=Wheel(10)):
    """Given a positive integer n, generate the list of primes <= n."""
    n += 1
    wS = wheel.size
    U = wheel.units
    nU = len(wheel.units)
    s = wheelSieveInner(n, wheel)
    return (wheel.primes[:bisect_left(wheel.primes, n)] +
            [p//nU*wS + U[p%nU] for p in xrange(bisect_left(U, 2 % wS)
             + 2//wS*nU, len(s)) if s[p]])

至于什么是轮子:嗯,你知道(除了 2),所有素数都是奇数,所以大多数筛子会漏掉所有偶数。同样,您可以走得更远一点,注意到所有素数(2 和 3 除外)都与 1 或 5 模 6 (== 2 * 3) 一致,因此您只需将这些数字的条目存储在筛子中即可. 下一步是注意所有素数(除了 2、3 和 5)都与 1、7、11、13、17、19、23、29(模 30)中的一个全等(这里 30 == 2*3 *5) 等等。

于 2010-05-24T14:57:59.887 回答
3

您可以使用位串进行的一项速度改进是在排除当前数字的倍数时使用“设置”方法。

所以要害部分变成

for i in range(3, limit, 2):
    if flags[i]:
        yield i
        if i <= sub_limit:
            flags.set(1, range(i*3, limit, i*2))    

在我的机器上,它的运行速度比原来的快 3 倍。

我的另一个想法是使用位串仅表示奇数。然后,您可以找到未设置的位,使用flags.findall([0])它返回一个生成器。不确定这是否会更快(找到位模式并不是很容易有效地完成)。

[完全披露:我编写了位串模块,所以我在这里有一些自豪感!]


作为比较,我还从 bitstring 方法中取出了胆量,以便它以相同的方式执行它,但没有范围检查等。我认为这为适用于十亿个元素的纯 Python 提供了一个合理的下限(没有改变算法,我认为这是作弊!)

def prime_pure(limit=1000000):
    yield 2
    flags = bytearray((limit + 7) // 8)
    sub_limit = int(limit**0.5)
    for i in xrange(3, limit, 2):
        byte, bit = divmod(i, 8)
        if not flags[byte] & (128 >> bit):
            yield i
            if i <= sub_limit:
                for j in xrange(i*3, limit, i*2):
                    byte, bit = divmod(j, 8)
                    flags[byte] |= (128 >> bit)

在我的机器上,一百万个元素的运行时间约为 0.62 秒,这意味着它大约是 bitarray 答案速度的四分之一。我认为这对于纯 Python 来说是相当合理的。

于 2010-05-24T13:53:25.380 回答
2

Python 的整数类型可以是任意大小,因此您不需要一个聪明的位串库来表示一组位,只需一个数字即可。

这是代码。它可以轻松处理 1,000,000 的限制,甚至可以处理 10,000,000 而不会抱怨太多:

def multiples_of(n, step, limit):
    bits = 1 << n
    old_bits = bits
    max = 1 << limit
    while old_bits < max:
        old_bits = bits
        bits += bits << step
        step *= 2
    return old_bits

def prime_numbers(limit=10000000):
    '''Prime number generator. Yields the series                                                                        
    2, 3, 5, 7, 11, 13, 17, 19, 23, 29 ...                                                                              
    using Sieve of Eratosthenes.                                                                                        
    '''
    yield 2
    sub_limit = int(limit**0.5)
    flags = ((1 << (limit - 2)) - 1) << 2
    # Step through all the odd numbers                                                                                  
    for i in xrange(3, limit, 2):
        if not (flags & (1 << i)):
            continue
        yield i
        # Exclude further multiples of the current prime number                                                         
        if i <= sub_limit:
            flags &= ~multiples_of(i * 3, i << 1, limit)

multiples_of功能避免了单独设置每个倍数的成本。它设置初始位,将其移动初始步骤 ( i << 1) 并将结果添加到自身。然后它将步长加倍,以便下一次移位复制两个位,依此类推,直到达到极限。这会将一个数字的所有倍数设置为 O(log(limit)) 时间的限制。

于 2010-05-24T13:44:01.440 回答
2

您可能想要查看的一个选项只是编译一个 C/C++ 模块,这样您就可以直接访问位旋转功能。AFAIK,您可以编写类似性质的东西,并且仅在完成用 C/C++ 执行的计算后返回结果。现在我正在输入此内容,您还可以查看 numpy,它确实将数组存储为本机 C 以提高速度。不过,我不知道这是否会比位串模块更快!

于 2010-05-24T14:29:45.837 回答
1

另一个非常愚蠢的选择,但如果你真的需要一个快速可用的大量素数列表,这可能会有所帮助。比如说,如果您需要它们作为回答项目 Euler 问题的工具(如果问题只是将您的代码优化为智力游戏,那就无关紧要了)。

使用任何缓慢的解决方案来生成列表并将其保存到 python 源文件中,说它只primes.py包含:

primes = [ a list of a million primes numbers here ]

现在要使用它们,您只需要以import primes磁盘 IO 的速度以最小的内存占用获得它们。复杂性是素数:-)

即使您使用优化不佳的解决方案来生成此列表,它也只会执行一次,并且没有太大关系。

您可能可以使用 pickle/unpickle 使其更快,但您明白了……

于 2010-05-27T00:36:23.097 回答
1

您可以使用分段的 Eratosthenes 筛。用于每个段的内存被重新用于下一个段。

于 2012-02-02T13:58:07.913 回答
1

这是一些 Python3 代码,它使用的内存比迄今为止发布的 bitarray/bitstring 解决方案少,大约是 Robert William Hanks 的 primesgen() 内存的 1/8,同时运行速度比 primesgen() 快(略快于 1,000,000,使用 37KB 内存, 在 1,000,000,000 时使用 34MB 以下,比 primesgen() 快 3 倍)。增加块大小(代码中的可变块)会使用更多内存,但会加速程序,达到一个限制——我选择了一个值,使其对内存的贡献低于筛子 n//30 字节的 10%。它不如Will Ness 的无限生成器(另请参见https://stackoverflow.com/a/19391111/5439078)内存效率高,因为它记录(并在最后以压缩形式返回)所有计算的素数。

只要平方根计算结果准确(如果 Python 使用 64 位双精度,大约为 2**51),这应该可以正常工作。但是你不应该使用这个程序来找到那么大的素数!

(我没有重新计算偏移量,只是从罗伯特威廉汉克斯的代码中获取它们。谢谢罗伯特!)

import numpy as np
def prime_30_gen(n):
    """ Input n, int-like.  Yields all primes < n as Python ints"""
    wheel = np.array([2,3,5])
    yield from wheel[wheel < n].tolist()
    powers = 1 << np.arange(8, dtype='u1')
    odds = np.array([1, 7, 11, 13, 17, 19, 23, 29], dtype='i8')
    offsets=np.array([[0,6,10,12,16,18,22,28],[6,24,16,12,4,0,22,10],
                      [0,6,20,12,26,18,2,8],  [24,6,4,18,16,0,28,10],
                      [6,24,26,12,14,0,2,20], [0,24,10,18,4,12,28,22],
                      [24,6,14,18,26,0,8,20], [0,24,20,18,14,12,8,2]],
                     dtype='i8')
    offsets = {d:f for d,f in zip(odds, offsets)}
    sieve = 255 * np.ones((n + 28) // 30, dtype='u1')
    if n % 30 > 1: sieve[-1] >>= 8 - odds.searchsorted(n % 30)
    sieve[0] &= ~1 
    for i in range((int((n - 1)**0.5) + 29) // 30):
        for odd in odds[(sieve[i] & powers).nonzero()[0]]:
            k = i * 30 + odd
            yield int(k)
            for clear_bit, off in zip(~powers, offsets[odd]): 
                sieve[(k * (k + off)) // 30 :: k] &= clear_bit
    chunk = min(1 + (n >> 13), 1<<15)
    for i in range(i + 1, len(sieve), chunk):
        a = (sieve[i : i + chunk, np.newaxis] & powers)
        a = np.flatnonzero(a.astype('bool')) + i * 8
        a = (a // 8 * 30 + odds[a & 7]).tolist()
        yield from a
    return sieve

旁注:如果您想要真正的速度并且将素数列表所需的 2GB 内存设置为 10**9,那么您应该使用 pyprimesieve(在https://pypi.python.org/上,使用 primesieve http:// primesieve.org/)。

于 2015-10-30T03:49:07.343 回答