有几个小的优化可以针对你的版本进行。通过颠倒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)
for i in range(3, limit, 2):
if flags[i]:
continue
yield i
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)
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
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)
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
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)
。使用以下代码,对于所有小于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)
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.'''
sieve_limit = gmpy2.isqrt(limit) + 1
limit += 1
bitmap = gmpy2.xmpz(3)
bitmap[4 : limit : 2] = -1
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 Honumpy.nonzero()
在这种情况下返回筛子数组中是质数的索引。你应该试一下。 - jfsgmpy
的版本比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