代码之家 › 专栏 › 技术社区 › Xavier Ho

加快python中的位串/位操作?

sieve-of-eratosthenes primes bit-manipulation optimization python-3.x

Xavier Ho · 技术社区 · 15 年前

我用 Sieve of Eratosthenes 以及python 3.1。代码以0.32秒的速度正常运行 ideone.com 生成最多1000000个质数。

# 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

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

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

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

您可以对这三行进行注释/取消注释,以查看我更改为使用位字符串的内容,如上面的代码片段所示。

我的问题是, 有没有一种方法可以加速我的程序,有没有位串?

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

再次编辑:

I've compiled a list of benchmarks on my machine.

11 回复 | 直到 6 年前

casevh 12 年前

您的版本有几个小的优化。通过颠倒“真”和“假”的角色,你可以改变” if flags[i] is False: “to” 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秒的运行时间。100000000的运行时间为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以接受 xmpz.bit_set(iterator) . 使用以下代码,所有小于100000000的素数的运行时间对于python 2.7为56秒,对于python 3.2为74秒。(如评论所述, xrange 比快 范围 )

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) . 使用下面的代码,所有小于100000000的素数的运行时间对于python 2.7和python 3.2都是大约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秒左右。(请参见gmpy2 2.0.0b1中的“修改”4。)

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

Xavier Ho 8 年前

好的,这里有一个(接近完成的)全面的基准测试,我今晚已经做了,看看哪个代码运行得最快。希望有人会发现这个列表有用。我省略了在我的机器上完成任何超过30秒的事情。

我要感谢每一个输入的人。我从你的努力中获得了很多见解,我希望你也有。

我的机器:AMDZM-86,2.40GHz双核,4GB内存。这是一台HP TouchSmart TX2笔记本电脑。请注意,虽然我可能链接到了一个Pastebin, 我在自己的机器上对以下所有内容进行了基准测试。

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

所有的基准测试都在python 2.6x86中进行。

返回质数n的列表,最多1000000个:( 使用蟒蛇发电机

Sebastian's numpy generator version (updated) - 121毫秒@

Mark's Sieve + Wheel - 154毫秒

Robert's version with slicing - 159毫秒

My improved version with slicing - 205毫秒

Numpy generator with enumerate - 249毫秒@

Mark's Basic Sieve - 317毫秒

casevh's improvement on my original solution - 343毫秒

My modified numpy generator solution - 407毫秒

My original method in the question - 409毫秒

Bitarray Solution - 414毫秒@

Pure Python with bytearray - 1394毫秒@

Scott's BitString solution - 6659 MS@

“@”表示此方法在我的机器设置。

另外,如果你不需要生成器,只需要整个列表马上:

numpy solution from RosettaCode - 32毫秒@

(numpy解决方案也能产生高达10亿,耗时61.6259秒。我怀疑记忆曾被交换过一次,因此是双倍的时间。)

Scott Griffiths 7 年前

好吧,这是我的第二个答案,但速度是关键,我认为我必须提到 bitarray 模块-即使它 bitstring “克星:”它非常适合这种情况,因为它不仅是C扩展(而且比纯Python希望的速度更快),而且还支持切片分配。但它还不能用于Python3。

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

10亿美元,运行良好,在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

James Camfield 6 年前

嗨,我正在寻找一个用python编写的代码来生成 10 ^ 9 尽可能快,这是因为记忆问题很困难。到目前为止,我想出了这个方法来生成素数。 10 ^ 6 和; 10 ^ 7 (在我的旧机器上分别计时0171s和1764s),这似乎比 “我改进的切片版本” 可能是因为我用了 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(标志[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

附言:如果您对如何加速这段代码有任何建议,从更改操作顺序到预计算,请发表评论。

Mark Dickinson Alexandru 15 年前

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

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

我有更快的版本,使用一个轮子,但它们更复杂。原始来源是 here .

好的,这是使用轮子的版本。 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),依此类推。

Scott Griffiths 15 年前

当排除当前数字的倍数时,可以使用“set”方法来提高使用位串的速度。

所以关键部分变成

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提供了一个合理的下限,它可以处理10亿个元素(不改变算法,我认为这是在欺骗!)

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秒的时间内运行一百万个元素,这意味着它大约是位数组应答速度的四分之一。我认为这对于纯Python来说是很合理的。

Marcelo Cantos 15 年前

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

这是密码。它可以轻松处理1000000个限额,甚至可以处理10000000个而不抱怨太多:

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))时间的限制。

Wayne Werner 15 年前

你可能想看的一个选项是编译C/C++模块,这样你就可以直接访问比特旋转特性。你可以写一些性质的东西,只返回在C/C++中完成计算的结果。现在我输入这个命令,您也可以看到numpy,它将数组存储为本机C以提高速度。不过,我不知道这是否比位串模块快!

kriss 15 年前

另一个非常愚蠢的选择,但是如果你真的需要一个大量的素数列表,那会很有帮助。比如说,如果你需要他们作为一个工具来回答ProjectEuler的问题(如果问题只是把你的代码优化成一个智力游戏,那是不相关的)。

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

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

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

即使您使用了一个不太优化的解决方案来生成这个列表,它也只会执行一次,而且并不重要。

你也许可以用泡菜/松饼来加快速度,但你明白了…

user448810 13 年前

你可以用埃拉托斯滕的分段筛。用于每个段的内存将在下一个段中重用。

Will Ness Derri Leahy 6 年前

这里有一些python3代码使用的内存比目前发布的bitarray/bitstring解决方案少,大约是robert william hanks的primesgen()的1/8内存,同时运行速度比primesgen()快(使用37kb内存时略快于1000000,使用34mb内存时略快于primesgen(),比primesgen()快3倍)。增加块大小(代码中的可变块)会使用更多的内存,但会加快程序的速度,直到达到一个限制——我选择了一个值,以便它对内存的贡献低于screen的n//30字节的10%。它的记忆效率不如 Will Ness's infinite generator (也见) https://stackoverflow.com/a/19391111/5439078 )因为它记录(并以压缩形式返回)所有计算出的素数。

只要平方根计算结果准确(如果python使用64位double,则约为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

旁注:如果你想要真正的速度,并且有2 GB的内存,那么你应该使用pyprimeshef(打开 https://pypi.python.org/ ,使用Prime谷粒筛 http://primesieve.org/ )