加快Python中的bitstring / bit操作？
我使用Eratosthenes Sieve和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
问题是，当我尝试产生高达10亿的数字时，内存不足。
flags = [False, False] + [True] * (limit - 2)
MemoryError
正如你可以想象的那样，分配10亿个布尔值（在Python中分配1个字节4或8个字节（见注释））实际上是不可行的，所以我研究了bitstring。 我想，每个标志使用1位会更有效率的记忆。 但是，该程序的性能急剧下降 - 运行时间为24秒，素数高达1,000,000。 这可能是由于bitstring的内部实现。
您可以评论/取消注释这三行，以查看我更改为使用BitString的内容，如上面的代码片段所示。
我的问题是， 有没有一种方法可以加快我的程序，无论有没有bittring？
编辑：请在发布之前亲自测试代码。 自然，我无法接受比现有代码慢的答案。
再次编辑：
我在我的机器上编译了一个基准测试列表。
您的版本有几个小的优化。 通过颠倒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来接受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))
好的，今天晚上我完成了一个（接近完整的）综合基准测试，以查看哪些代码运行速度最快。 希望有人会发现这个列表有用。 我省略了需要30秒以上完成的任何事情。
我想感谢大家提出意见。 我从你的努力中获得了很多见解，我也希望你也有。
我的机器：AMD ZM-86,2.40 Ghz双核，带有4GB内存。 这是惠普Touchsmart Tx2笔记本电脑。 请注意，虽然我可能已经链接到一个pastebin，但我在我自己的机器上测试了以下所有内容。
一旦我能够构建它，我将添加gmpy2基准测试。
所有的基准测试都在Python 2.6 x86中进行测试
返回素数列表n（高达1,000,000：）（使用Python生成器）
塞巴斯蒂安的numpy发电机版本（更新） - 121毫秒@
马克的筛+车轮 - 154毫秒
罗伯特的版本与切片 - 159毫秒
我的改进版与切片 - 205毫秒
具有枚举的Numpy生成器 - 249 ms @
马克的基本筛 - 317毫秒
casevh改进了我的原始解决方案 - 343毫秒
我修改后的numpy发生器解决方案 - 407毫秒
我在问题中的原始方法 - 409毫秒
Bitarray解决方案 - 414 ms @
纯Python与bytearray - 1394毫秒@
Scott的BitString解决方案 - 6659 ms @
'@'表示这种方法能够在我的机器设置上产生高达n <1,000,000,000的数据。
另外，如果你不需要生成器，只需要一次列出整个列表：
来自RosettaCode的numpy解决方案 - 32 ms @
（numpy的解决方案也能产生高达10亿次，这花费了61.6259秒，我怀疑内存换了一次，因此是双倍时间。）
好吧，这是我的第二个答案，但是由于速度是本质的，我认为我必须提到bitarray模块 - 尽管它是位串的克星:)。 它非常适合这种情况，因为它不仅是一个C扩展（比纯Python有希望存在的速度更快），而且还支持片分配。 尽管如此，它还不适用于Python 3。
我甚至没有试图优化这个，我只是重写了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 xrange(3, limit, 2):
if not flags[i]:
yield i
if i <= sub_limit:
flags[3*i:limit:i*2] = True