加快 Python 中的位串/位操作？
Posted
技术标签:
【中文标题】加快 Python 中的位串/位操作？【英文标题】：Speed up bitstring/bit operations in Python? 【发布时间】：2011-02-23 05:57:03 【问题描述】：我使用Sieve of Eratosthenes 和 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 亿个布尔值（1 字节在 Python 中每个 4 或 8 个字节（见注释））确实不可行，所以我查看了 bitstring。我想，为每个标志使用 1 位会更节省内存。然而，该程序的性能急剧下降 - 运行时间为 24 秒，质数高达 1,000,000。这可能是由于位串的内部实现。
您可以注释/取消注释这三行以查看我更改为使用 BitString 的内容，如上面的代码 sn-p。
我的问题是，有没有办法加快我的程序，有或没有位串？
编辑：请在发布前自行测试代码。我自然不能接受运行速度比我现有代码慢的答案。
再次编辑：
I've compiled a list of benchmarks on my machine.
【问题讨论】：
Nitpick：您正在查看每个布尔值 4 或 8 个字节（取决于您是在 32 位还是 64 位系统上），而不是 1：在内部，列表 @987654328 @ 只是一个 (PyObject *) 指针的 C 数组。 感谢指正。:]
你可以在 Python 2.x 中使用 numpy rosettacode.org/wiki/Sieve_of_Eratosthenes#Using_numpy 它更快（~20 倍）。
使用上述链接中的 primes_upto2()、primes_upto3() 生成最多限制 1e9 需要大约 28 秒。 C++ 版本 ~14 秒，C 版本 ~13 秒。
@wallacaloo：在 Python 3.x 中，范围是惰性的。
【参考方案1】：
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)) 时间的限制。
【讨论】：
@Marcelo：是的，我知道这一点。我还可以编写一些执行位操作的函数，看看它是否更快。这是我目前正在做的事情。编辑：问题是，即使执行2 << 1000000 也需要超过 10 秒。
@Xavier：不，它没有。打印该操作的结果可能需要很长时间。请改用x = 2 << 1000000。
@kaizer.se：嗯，好点子。我去看看我能用它做什么。
这似乎不太可能奏效，除非你能找到某种方法将一大堆位集放在一起：long 上的每个位操作都会生成一个全新的 long，O 也是(n) 操作。
@Mark：是的，我刚试过。 1.3 秒生成最多 100,000，当我尝试 1,000,000 时需要 20 多秒（实际上，我键入时它还没有完成）。这是通过位移不可行的，我需要直接访问来更改位。【参考方案2】：
使用位串可以提高速度的一个方法是在排除当前数字的倍数时使用“设置”方法。
所以要害部分变成了
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 来说是相当合理的。
【讨论】：
@Scott：很酷，很高兴这里有位字符串的作者！我也试试，很快就会告诉你结果。是的，我使用的是 2.0.0 beta 1，因为我一个小时前才下载了这个库。 @Scott：刚刚做了一个测试。 11.2 秒与您的修复。还是有点慢。还有什么想法吗？ 在上面加上更多的想法。我猜这应该会让你的时间减少到大约 7 秒。 @Scott：谢谢。但我的原始代码运行时间为 0.32 秒。见：ideone.com/wCJm5。但感谢您的想法和意见！我会关注这个话题一段时间。 我认为挑战在于最快的筛选代码产生多达 10 亿个素数，而不是 100 万个。位串代码可以工作十亿，而您的原始代码不能，所以我认为位串到目前为止是赢了！顺便说一句，使用 Python 2 我将百万案例缩短到 4.5 秒。【参考方案3】：您可能想要查看的一个选项只是编译一个 C/C++ 模块，这样您就可以直接访问位旋转功能。 AFAIK，您可以编写具有这种性质的东西，并且仅在完成用 C/C++ 执行的计算后才返回结果。现在我正在输入此内容，您还可以查看 numpy，它确实将数组存储为本机 C 以提高速度。不过，我不知道这是否会比位串模块更快！
【讨论】：
谢谢，韦恩。你是对的，它是一种选择——只是不是一个简单的选择。当然，当你写一个时，我会很高兴。;]
嘿。如果您已经了解 C/C++，这实际上并没有那么糟糕——最大的痛苦是弄清楚如何告诉您的脚本正确的事情，让 scons 来处理构建过程。然后你只需要处理大约 125 MB 的比特（10 亿比特/8 字节 == 1.25 亿字节）。
是的。我知道 C++，但我知道 Python 是用 C 编写的，而且我自己还没有用 C/C++ 编写过 Python 模块。所以这还有些遥远。不过没关系，我们一如既往地在 SO 上得到了一些出色的答案。 :]【参考方案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
我有更快的版本，使用***，但它们要复杂得多。原始来源是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) 等等。
【讨论】：
愿意解释一下吗？是不是类似于……阿金筛子？ @Mark：0.28 秒！你是我们第一个进入决赛的人！ =Dideone.com/yIENn @Mark：当然，MemoryError@1,000,000,000。 =/ideone.com/YeBOR。我现在对你的更快版本很好奇。
谢谢！我将阅读代码并尝试理解它。稍后会报告我的状态。
不，阿特金的 Sieve 引入了一个全新的想法，即使用二次形式；我认为它只是渐近地比一个基本的eratosthenes + wheel 筛子快，并且它变得更快的点可能是> 1000000（当然取决于实现）。使用***的想法已经有很长一段时间了。我添加了一个链接到我第一次发布这个的地方；那里有一个使用***的实现。【参考方案5】：
好的，这是我的第二个答案，但由于速度至关重要，我认为我不得不提及 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
【讨论】：
哇什么，位数组！正是我需要的？ XD。我今天下班后试试看。 是的，我之前也遇到过同样的问题，并打算建议使用 bitarray。我以前没有听说过比特串，但我会检查一下。在我发现 bitarray 之前我一直在使用 BitVector（我发现它不是很优化），但是很高兴知道另一个二进制数组模块要检查。 只是想让你知道，在我的机器上，运行 n 【参考方案6】：您的版本有几个小的优化。通过颠倒 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 秒。 （如 cmets 中所述，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))
【讨论】：
你用什么设置来安装gmpy。在我的机器上 limit=1e9 需要 500 秒（相比之下，rosettacode.org/wiki/Sieve_of_Eratosthenes#Using_numpy 的 primes_upto2() 需要 30 秒）。我刚刚检查了代码并运行python setup.py install
@casevh：感谢 gmpy 代码。今天下班后我会试一试 - 减少 2/3 的运行时间给我留下了深刻的印象。
这应该是您安装它所需要做的所有事情。我正在使用 64 位 Ubuntu 10.04、2.2Ghz Core2 Duo、GMP 5.01 和 Python 3.1。使用 Ubuntu 版本的 numpy，primes_upto2() 在我的计算机上需要 50 秒，所以有些奇怪。
@casevh：我注意到prime_numbers() 的第二个变体导致我的机器发生交换。所以它可能会破坏基准。 64 位 Ubuntu 10.04、2.6GHz i7、GMP 4.3.2 pgmpy 2.0.0a 和 python 2.6.4。
由于您使用的是 Python 2.x，请将范围更改为 xrange。那应该可以解决极端的内存使用问题。您可以尝试从源代码编译 GMP 5.01。 IIRC，它为较新的英特尔处理器改进了代码。【参考方案7】：
好的，这是我今晚完成的（接近完整的）综合基准测试，以查看哪些代码运行速度最快。希望有人会发现这个列表很有用。我省略了在我的机器上需要超过 30 秒才能完成的任何内容。
我要感谢所有提出意见的人。我从你的努力中获得了很多见解，我希望你也一样。
我的机器：AMD ZM-86，2.40 Ghz 双核，4GB RAM。这是一台 HP Touchsmart Tx2 笔记本电脑。请注意，虽然我可能已链接到 pastebin，但我在自己的机器上对以下所有内容进行了基准测试。
一旦我能够构建它，我将添加 gmpy2 基准。
所有基准测试均在 Python 2.6 x86 中测试
返回一个质数列表 n 最多 1,000,000：（使用 Python 发电机）
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 女士@
'@' 表示此方法最多可以生成 n
此外，如果您不需要 生成器，只想要整个列表 一次：
numpy solution from RosettaCode - 32 毫秒@
（numpy 的解决方案也能够生成多达 10 亿个，耗时 61.6259 秒。我怀疑内存被交换了一次，因此是双倍的时间。）
【讨论】：
@Xavier Ho：你的 numpy 版本不正确：步骤应该是n，而不是n*n，即isprime[n*n::n]=0。如果您想要生成器版本，则无需致电 numpy.nonzero()：primes[:2]=0; return (i for i, p in enumerate(primes) if p)
注意：numpy 生成器解决方案比非生成器版本慢 3 倍（100 秒对 1e9 的 30 秒）。
@J.F Sebastian：不错的收获。以为我修好了！谢谢。
@J.F.塞巴斯蒂安：有趣。在我的机器上，它慢了 6 倍以上。
@Xavier Ho：这里是 numpy 生成器版本ideone.com/mxSdH（1e9 为 55 秒（相比之下，numpy 非生成器版本为 30 秒，@casevh 基于 xmpz.bitset() 的新版本为 45 秒）版本）【参考方案8】：
相关问题：Fastest way to list all primes below N in python。
嗨，我也在寻找 Python 中的代码以尽可能快地生成高达 10^9 的素数，由于内存问题，这很困难。到目前为止，我想出这个来生成最高 10^6 和 10^7 的素数（在我的旧机器上分别计时 0,171s 和 1,764s），这似乎比上一篇文章中的 “我的切片改进版” 稍微快一点（至少在我的电脑上），可能是因为我使用 n//i-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：如果您对如何加快此代码的速度有任何建议，从更改操作顺序到预先计算任何内容，请发表评论。
【讨论】：
嗯，它会更快，因为您使用的是列表理解而不是生成器。不错，有时间我会加基准的。 只是一个想法，您能解释一下您的m = n // i-i 与我的flags[i2::i<<1] 有何相似之处吗？由于我忽略了遍历所有 2 的倍数，因此我也在切片语法的“步骤”中避免了它。你这样做的方式仍然是迭代自身的每一个倍数。
对不起，我是编程新手，我什至不知道 i）中的倍数 oi 的数量是 i (1i,2i,3i,...ii) 所以在 [ii::i] 我们有 i 的倍数在 range(1,n) - i 的倍数在 range(1, ii) +1 (+one 因为我们也想计算 ii) 因此公式 len(sieve[i*i::i]) 等于 n//i-i +1。
在我的代码中，我完全忽略了 2 的倍数，并且我没有标记出我的筛子，因为我依靠 range 函数的步骤为 2 来筛选和形成素数列表。（我仅将 [2] 添加到列表末尾）
另一方面，如果筛子的初始化是这样的 [False,True,True] * ((n+1)//3) ，也可以完全跳过素数 3 见我的第二个例子 primes2()，它有点快。请确保输入小于 3 的倍数。在我的机器中，差异是如此之小，以至于我不关心更好的代码。【参考方案9】：
另一个非常愚蠢的选择，但如果你真的需要一个快速可用的大量素数列表，这可能会有所帮助。比如说，如果您需要它们作为解决项目 Euler 问题的工具（如果问题只是将您的代码优化为智力游戏，那就无关紧要了）。
使用任何缓慢的解决方案来生成列表并将其保存到 python 源文件中，primes.py 说，它只包含：
primes = [ a list of a million primes numbers here ]
现在要使用它们，您只需执行import primes 即可以磁盘 IO 的速度以最少的内存占用获得它们。复杂性是素数的数量:-)
即使您使用优化不佳的解决方案来生成此列表，它也只会执行一次，并且没有太大关系。
您可能可以使用 pickle/unpickle 使其更快，但您明白了...
【讨论】：
【参考方案10】：您可以使用分段的埃拉托色尼筛。每个段使用的内存被下一个段重用。
【讨论】：
“分段”是指某个数字范围的内存块，一旦用完，您就在同一内存块上创建下一个数字范围？ 是的。谷歌的“埃拉托色尼分段筛”。 s/可以/应该。 :)【参考方案11】：这是一些 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's infinite generator（另见https://***.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 的 RAM 来将素数列表设置为 10**9，那么您应该使用 pyprimesieve（在 https://pypi.python.org/ 上，使用 primesieve http://primesieve.org/）。
【讨论】：
以上是关于加快 Python 中的位串/位操作？的主要内容，如果未能解决你的问题，请参考以下文章
python 将一个转换为一个位数组字符串，看看你在python中的位操作中做了什么。