如何在Python中加速比特串/比特操作?

44

我使用埃拉托斯特尼筛法和Python 3.1编写了一个质数生成器。该代码在ideone.com上运行,优雅地正确生成小于1,000,000的质数,所需时间为0.32秒。

# 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

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

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

正如您所想象的那样,分配10亿个布尔值(在Python中每个值占用4或8个字节)是不可行的,因此我研究了bitstring。我想,为每个标志使用1位将更加节省内存。然而,程序的性能急剧下降-运行时间为24秒,适用于小于100万的质数。这可能是由于bitstring的内部实现造成的。
您可以取消/注释三行以查看我更改使用BitString的内容,如上面的代码片段。
我的问题是,是否有一种方法可以加速我的程序,无论是否使用bitstring? 编辑:请在发布之前自行测试代码。我不能接受比我现有代码运行得更慢的答案,这是自然的。 再次编辑: 我已经在我的机器上编译了一份基准列表。

使用上述链接中的primes_upto2()和primes_upto3()函数生成1e9的极限需要约28秒。C++版本大约需要14秒,C版本则需要13秒。 - jfs
@J.F. Sebastian:不错的链接。今天下班后我会尝试一下,但要注意的是它没有使用生成器,并且仅返回一个由零和一组成的数组。并不完整。:] - Xavier Ho
Xavier Ho:numpy.nonzero()在这种情况下返回筛子数组中质数的索引。你应该试一下。 - jfs
casevh基于gmpy的版本比iprimes_upto()稍微快一点,但比基于numpy的版本慢得多。对于1e6:iprimes_upto: 0.39,primes_upto2: 0.02,primes_upto3: 0.02,prime_numbers: 0.29。对于1e9,prime_numbers太慢了(但没有内存错误)。primes_upto2: 28.96,primes_upto3: 27.55,prime_numbers: 509.74(以秒为单位的时间)。 - jfs
5
在Python 3.x中,range是懒惰的。 - jfs
显示剩余11条评论
11个回答

35

有几个小的优化可以针对你的版本进行。通过颠倒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秒。限制为10亿时,运行时间为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)。使用以下代码,对于所有小于1,000,000,000的质数,Python 2.7的运行时间为56秒,Python 3.2的运行时间为74秒。(如评论中所述,xrangerange更快。)

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)。使用以下代码,对于所有小于1,000,000,000的质数,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))

@casevh:Python 3.2a0 稍微慢了一些:prime_numbers: 173.18 prime_numbers2: 110.99 - jfs
@J.F. Sebastian:我的结果与您在Py3与Py2方面的观点一致。即使是Python 2.6中的xrange()也比Python 3中的range更快。我不知道为什么。 - Xavier Ho
@J.F. Sebastian:你在Python 3.2a0上的性能结果很有趣。我通常使用3.2获得相同或稍快的时间。为确保配置选项相同,我编译了每个用于测试的Python版本。你在编译3.2时是否使用了--with-computed-gotos?如果我没记错,它默认是禁用的。 - casevh
哇,世界变得越来越好了。干得好,casevh。我下周应该会安装gmp在我的电脑上,然后我会告诉你基准测试结果。 - Xavier Ho
第三次编辑的时间:Python 2.6:24秒,3.2:33秒,超过了原始的numpy版本http://rosettacode.org/wiki/Sieve_of_Eratosthenes#Using_numpy。整数的映射接口对我来说可能太神奇了。我将开启一个问题以解决遇到的问题。 - jfs
显示剩余23条评论

13

好的,这是我的第二个回答,但是由于速度至关重要,我认为我必须提到 bitarray 模块 - 尽管它是 bitstring 的劲敌 :). 它非常适合这种情况,因为它不仅是 C 扩展(快于纯 Python),而且还支持切片赋值。

我甚至没有尝试优化它,只是重写了 bitstring 版本。在我的电脑上,对于小于一百万的质数,我得到了 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 range(3, limit, 2):
        if not flags[i]:
            yield i
            if i <= sub_limit:
                flags[3*i:limit:i*2] = True

哇,位数组!正是我需要的东西?XD。今天下班后我会试一试。 - Xavier Ho
1
是的,我之前遇到过同样的问题,打算建议使用bitarray。我之前没有听说过bitstring,但我会去看看它。在发现bitarray之前,我一直在使用BitVector(我发现它并不是很优化),但知道还有另一个二进制数组模块可以尝试真是太好了。 - Justin Peel
只是想让你知道,在我的机器上,当 n < 1,000,000 时运行所需的时间为0.45秒,因此执行十亿需要大约450秒。我还没有尝试过,但是这可以让你对我的机器速度相对于我的0.21秒版本有些了解。也许我可以使用bitarray来生成需要超过1亿的东西,呵呵。 - Xavier Ho

8

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

我要感谢所有提供意见的人。我从你们的努力中获得了很多见解,希望你也是如此。

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

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

所有基准测试都在Python 2.6 x86中进行测试

使用 Python 生成器返回小于1,000,000的质数列表:

Sebastian的numpy生成器版本(已更新) - 121毫秒 @

Mark的筛法+轮子 - 154毫秒

Robert的切片版本 - 159毫秒

我改进的切片版本 - 205毫秒

带有枚举的Numpy生成器 - 249毫秒 @

Mark的基本筛法 - 317毫秒

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

我修改过的Numpy生成器解决方案 - 407毫秒

我的问题中的原始方法 - 409毫秒

Bitarray解决方案 - 414毫秒 @

使用bytearray的纯Python - 1394毫秒 @

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

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

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

RosettaCode中的numpy解决方案 - 32毫秒 @

(Numpy解决方案也能够生成高达10亿,用时61.6259秒。我怀疑内存被交换了一次,因此时间加倍。)


@J.F. Sebastian:有趣。在我的机器上,它比你的慢了6倍以上。 - Xavier Ho
ambi_sieve()比RossettaCode上的numpy解决方案更快。https://dev59.com/snI95IYBdhLWcg3w-DH0#2068548 - jfs
primesfrom2to()在n=1e9时需要10秒钟。这是我所知道的Python解决方案中最好的结果。https://dev59.com/snI95IYBdhLWcg3w-DH0#3035188 - jfs
@Xavier Ho:我添加了我最后一次尝试的numpy生成器,如果你有时间,可以使用timeit测试一下。(我知道这是一个维基,但在其他机器上进行一些基准测试会很好)。 - Robert William Hanks
我并不建议重新唤起这个古老的线程,但是如果没有在几个规模点比较运行时间,并获取代码的经验增长率,就无法进行适当的速度比较(ideal情况下),并将时间和大小的图形绘制在一个对数-对数图上,例如这里这里所示。 - Will Ness
显示剩余10条评论

6

相关问题:Python中列出所有小于N的质数最快的方法

你好,我也在寻找一种能够尽可能快地生成小于10^9的质数的Python代码,但由于内存问题这很困难。目前为止,我已经写出了以下代码来生成小于10^610^7的质数(在我的旧机器上分别需要0.171秒和1.764秒),它似乎比之前帖子中的"我的改进版本"稍微快一些(至少在我的电脑上是这样),可能是因为我使用了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*i2*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:如果您对如何加速此代码有任何建议,从更改操作顺序到预计算任何内容,都请留言。


好的,这样会更快,因为你使用了列表推导而不是生成器。很棒,我有时间会添加基准测试。 - Xavier Ho
对不起,我刚开始学编程,现在甚至不知道 << 是什么意思。而且我也不确定我的代码是否比你的快,或者为什么快,我猜可能是你调用了 len()。也许这可以帮到你,如果离题了请原谅。嗯,数学告诉我们,在范围(1,n)内,i的倍数的个数等于 n//i(整数除法),在范围(1,ii)内,i的倍数的个数为 i (1i,2i,3i,...ii) ,所以在[ii::i]中,我们有范围(1,n)内的i的倍数 - 范围(1,ii)内的i的倍数 +1(加一是因为我们还要计算 i*i)因此公式len(sieve[i**i::i])等于n//i-i+1。 - Robert William Hanks
顺带一提,如果筛子的初始化像我的第二个例子primes2()那样是这样的[False,True,True] * ((n+1)//3),人们也可以完全跳过3这个质数,它会更快一些。请确保输入比3的倍数少1。在我的机器上,差别非常小,以至于我不需要更好的代码。 - Robert William Hanks
primesgen2() 相较于 primesgen1() 有微小的改进:对于 n=1e6,时间分别为 33 毫秒和 37 毫秒。 http://ideone.com/xcoii - jfs
当n=1亿时,结果证实了“内存假说”:primesgen1()primesgen2()的运行时间分别为35秒和27秒(时间比例与内存使用比例相似)。而primesfrom2to3()则只需要9.9秒(作为对比)。 - jfs
显示剩余9条评论

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)模6余1或5(== 2*3),因此你只需要在筛子中存储那些数字的条目。再往上一步,你会发现所有的质数(除2、3和5外)模30余1、7、11、13、17、19、23或29(这里30 == 2*3*5)。

能否解释一下“whee”的含义?它是否类似于...Akin筛法? - Xavier Ho
@Mark:0.28秒!你是我们晋级决赛的第一名!=D http://ideone.com/yIENn - Xavier Ho
@Mark:自然而然,MemoryError @ 1,000,000,000。=/ http://ideone.com/YeBOR。我现在很好奇你的更快版本。 - Xavier Ho
谢谢! 我会阅读代码并尝试理解它。稍后会报告我的状态。 - Xavier Ho
不,Atkin筛法引入了一种根本新颖的想法,即使用二次形式;我认为它只在*渐进意义下比Eratosthenes筛法+Wheel快,并且它变得更快的点可能会大于1000000(当然取决于实现)。使用轮子的想法已经存在了一段时间。我已经添加了一个链接到我最初发布这个内容的地方;那里有一个使用轮子的实现。 - Mark Dickinson
显示剩余2条评论

3
使用bitstring可以提高速度的一种方法是,在排除当前数字的倍数时使用“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模块,所以我在这里有些自豪!]
作为比较,我还取出了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来说是相当合理的。


@Scott:噢,太酷了,很高兴在这里有bitstring的作者!我也会尝试一下,很快会告诉你结果的。是的,我正在使用2.0.0 beta 1版本,因为我只是在一个小时之前下载了这个库。 - Xavier Ho
@Scott:刚刚做了一个测试。使用你的修复后,需要11.2秒。仍然有点慢。你还有什么更好的想法吗? - Xavier Ho
使用筛法算质数对我来说感觉像是作弊 - 它不再是同样的算法了。通过剥离一些内部位串验证代码,我将计算质数到一百万的时间缩短到了2.3秒。等它计算十亿的时候我会告诉你需要多长时间 :) - Scott Griffiths
@Scott:感谢你提供的并排比较。是的,你做得很好。:3 - Xavier Ho
prime_pure()的速度比我比较过的最慢的函数慢两倍 http://ideone.com/F0RBo(比较在文件末尾) - jfs
显示剩余6条评论

2

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

以下是代码。它可以轻松处理100万的限制,甚至可以处理1000万而不会抱怨太多:

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))时间内设置了所有小于限制的数字的倍数。


2
@Xavier:不是的。可能是打印该操作结果需要那么长时间。尝试使用x = 2 << 1000000 - u0b34a0f6ae
随着极限的增长,减速是不可避免的。只是出于好奇:bitset版本是更快还是更慢?你见过最快的100万次运行时间是多少? - Marcelo Cantos
@Marcelo:你指的是哪个bitset版本?如果你是指Scott的回答,那么在该版本中,当n < 1,000,000时,它需要11秒。尽管如此,这仍然比我的代码慢10倍。 - Xavier Ho
@Marcelo:我在这里找到了一段C#代码:http://ideone.com/1GMyI || 它声称可以在190毫秒内运行,并找到100万以下所有质数的总和。虽然我还没有测试过,但我相信它 - 这正是我的算法。 - Xavier Ho
当然,编译型语言比Python更快地解决这个问题。如果有一个用C编写的位串库(最好支持设置一条位串),那么你就可以兼得两种优势。 - Marcelo Cantos
显示剩余12条评论

2

你可能想要考虑的一种选择是编译一个C/C++模块,以便直接访问位操作功能。据我所知,你可以编写这样的代码,只需在C/C++中完成计算后返回结果即可。现在我正在打字时,你还可以看看numpy,它将数组存储为本地C,以提高速度。不过,我不知道这是否比bitstring模块更快!


谢谢,韦恩。你说得对,这是一个选项 - 只是不是很容易的选项。当然,如果你写一个,我会很高兴的。;] - Xavier Ho
嘿,如果你已经了解C/C++,那么实际上并不难 - 最大的痛点是弄清楚如何告诉你的脚本正确的事情,以便scons来处理构建过程。然后你只需要处理大约125 MB的位(10亿位/8字节== 1.25亿字节)。 - Wayne Werner
没错,我懂C++,但是我知道Python是用C写的,而且我自己还没有用C/C++编写过Python模块。所以这有点遥远。不过没关系,我们在SO上总能得到一些很棒的答案,一如既往。:] - Xavier Ho

1
这里有一些Python3代码,比到目前为止发布的bitarray/bitstring解决方案使用更少的内存,比Robert William Hanks的primesgen()占用的约1/8的内存更快运行,而且在1,000,000时仍然比primesgen()快(使用37KB的内存),在1,000,000,000时比primesgen()快3倍,但是chunk变量增加会增加程序的内存使用,不过有一个限制——我选择了一个值,以便其对筛法的n // 30字节的内存贡献不超过10%。它不像Will Ness的无限生成器(也参见https://dev59.com/K3E95IYBdhLWcg3wp_yn#19391111)那样内存效率高,因为它记录(并在压缩形式下返回)所有已计算的素数。
只要平方根计算准确(如果Python使用64位双精度,则大约为2 ** 51),这应该能正确运行。但是,您不应该使用此程序来查找如此大的素数!

我没有重新计算偏移量,而是使用了Robert William Hanks的代码。谢谢Robert!

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的RAM以生成10**9以下的质数列表,那么您应该使用pyprimesieve(在https://pypi.python.org/上,使用primesieve http://primesieve.org/)。


1

另一个非常愚蠢的选项,但如果您确实需要快速获取大量质数列表,则可能会有所帮助。比如说,如果您需要它们作为回答项目欧拉问题的工具(如果问题只是将代码优化为头脑游戏,则无关紧要)。

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

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

现在,只需执行import primes,即可使用它们,并以磁盘IO的速度和最小的内存开销获得质数。其复杂度就是质数个数 :-)

即使您使用了一个效率较低的解决方案来生成这个列表,它也只会完成一次,并且并不重要。

您可以尝试使用pickle/unpickle来进一步提高速度,但是您已经掌握了主要思路...


网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接