Как реализовать эффективный бесконечный генератор простых чисел в Python?
Это не домашняя работа, мне просто любопытно.
БЕСКОНЕЧНОЕ ключевое слово здесь.
Я хочу использовать это как для p в простых числах (). Я считаю, что это встроенная функция в Haskell.
Таким образом, ответ не может быть таким наивным, как "Просто сделай сито".
Прежде всего, вы не знаете, сколько последовательных простых чисел будет использовано. Ну, предположим, вы могли бы придумать 100 из них одновременно. Будете ли вы использовать тот же подход Сита, а также формулу частоты простых чисел?
Я предпочитаю не одновременный подход.
Спасибо, что читаете (и пишете;))!
15 ответов
"Если бы я видел дальше..."
erat2 Функция из поваренной книги может быть дополнительно ускорена (примерно на 20-25%):
erat2a
import itertools as it
def erat2a( ):
D = { }
yield 2
for q in it.islice(it.count(3), 0, None, 2):
p = D.pop(q, None)
if p is None:
D[q*q] = q
yield q
else:
# old code here:
# x = p + q
# while x in D or not (x&1):
# x += p
# changed into:
x = q + 2*p
while x in D:
x += 2*p
D[x] = p
not (x&1) проверка подтверждает, что x странно Тем не менее, так как оба q а также p странные, добавив 2*p половину шагов избегают вместе с тестом на странность.
erat3
Если кто-то не против немного лишней фантазии, erat2 может быть увеличен на 35-40% со следующими изменениями (примечание: требуется Python 2.7+ или Python 3+ из-за itertools.compress функция):
import itertools as it
def erat3( ):
D = { 9: 3, 25: 5 }
yield 2
yield 3
yield 5
MASK= 1, 0, 1, 1, 0, 1, 1, 0, 1, 0, 0, 1, 1, 0, 0,
MODULOS= frozenset( (1, 7, 11, 13, 17, 19, 23, 29) )
for q in it.compress(
it.islice(it.count(7), 0, None, 2),
it.cycle(MASK)):
p = D.pop(q, None)
if p is None:
D[q*q] = q
yield q
else:
x = q + 2*p
while x in D or (x%30) not in MODULOS:
x += 2*p
D[x] = p
erat3 Функция использует тот факт, что все простые числа (кроме 2, 3, 5) по модулю 30 приводят только к восьми числам: те, которые включены в MODULOS frozenset. Таким образом, после получения начальных трех простых чисел, мы начинаем с 7 и работаем только с кандидатами.
Фильтрация кандидатов использует itertools.compress функция; "магия" находится в MASK последовательность; MASK имеет 15 элементов (есть 15 нечетных чисел в каждых 30 числах, выбранных itertools.islice функция) с 1 для каждого возможного кандидата, начиная с 7. Цикл повторяется, как указано в itertools.cycle функция.
Введение фильтрации кандидатов требует еще одной модификации: or (x%30) not in MODULOS проверять. erat2 алгоритм обработал все нечетные числа; теперь, когда erat3 алгоритм обрабатывает только кандидатов R30, мы должны убедиться, что все D.keys() могут быть только такие "ложные" кандидаты.
Ориентиры
Результаты
На сервере Atom 330 Ubuntu 9.10 версии 2.6.4 и 3.1.1+:
$ testit
up to 8192
==== python2 erat2 ====
100 loops, best of 3: 18.6 msec per loop
==== python2 erat2a ====
100 loops, best of 3: 14.5 msec per loop
==== python2 erat3 ====
Traceback (most recent call last):
…
AttributeError: 'module' object has no attribute 'compress'
==== python3 erat2 ====
100 loops, best of 3: 19.2 msec per loop
==== python3 erat2a ====
100 loops, best of 3: 14.1 msec per loop
==== python3 erat3 ====
100 loops, best of 3: 11.7 msec per loop
На домашнем сервере AMD Geode LX Gentoo Python 2.6.5 и 3.1.2:
$ testit
up to 8192
==== python2 erat2 ====
10 loops, best of 3: 104 msec per loop
==== python2 erat2a ====
10 loops, best of 3: 81 msec per loop
==== python2 erat3 ====
Traceback (most recent call last):
…
AttributeError: 'module' object has no attribute 'compress'
==== python3 erat2 ====
10 loops, best of 3: 116 msec per loop
==== python3 erat2a ====
10 loops, best of 3: 82 msec per loop
==== python3 erat3 ====
10 loops, best of 3: 66 msec per loop
Контрольный код
primegen.py модуль содержит erat2, erat2a а также erat3 функции. Здесь следует сценарий тестирования:
#!/bin/sh
max_num=${1:-8192}
echo up to $max_num
for python_version in python2 python3
do
for function in erat2 erat2a erat3
do
echo "==== $python_version $function ===="
$python_version -O -m timeit -c \
-s "import itertools as it, functools as ft, operator as op, primegen; cmp= ft.partial(op.ge, $max_num)" \
"next(it.dropwhile(cmp, primegen.$function()))"
done
done
Поскольку ОП требует эффективной реализации, вот существенное улучшение кода 2002 активного состояния Дэвида Эппштейна / Алекса Мартелли (см. Здесь в его ответе): не записывайте информацию простого числа в словаре, пока его квадрат не будет виден среди кандидаты. Уменьшает сложность пространства до O (sqrt (n)) вместо O (n), для n произведенных простых чисел ( π (sqrt (n log n)) ~ 2 sqrt (n log n) / log (n log n) ~ 2 sqrt (n / log n)). Следовательно, временная сложность также улучшается, то есть она работает быстрее.
Создает "скользящее сито" в качестве словаря текущих кратных значений каждого базового простого числа (т. Е. Ниже квадрата текущей производственной точки) вместе с их значениями шага:
from itertools import count
# ideone.com/aVndFM
def postponed_sieve(): # postponed sieve, by Will Ness
yield 2; yield 3; yield 5; yield 7; # original code David Eppstein,
sieve = {} # Alex Martelli, ActiveState Recipe 2002
ps = postponed_sieve() # a separate base Primes Supply:
p = next(ps) and next(ps) # (3) a Prime to add to dict
q = p*p # (9) its sQuare
for c in count(9,2): # the Candidate
if c in sieve: # c's a multiple of some base prime
s = sieve.pop(c) # i.e. a composite ; or
elif c < q:
yield c # a prime
continue
else: # (c==q): # or the next base prime's square:
s=count(q+2*p,2*p) # (9+6, by 6 : 15,21,27,33,...)
p=next(ps) # (5)
q=p*p # (25)
for m in s: # the next multiple
if m not in sieve: # no duplicates
break
sieve[m] = s # original test entry: ideone.com/WFv4f
(старый, оригинальный код здесь был отредактирован, чтобы включить изменения, как видно из ответа Tim Peters ниже). см. также это для соответствующего обсуждения.
Подобный код на основе колес 2-3-5-7 работает примерно в 2,15 раза быстрее (что очень близко к теоретическому улучшению 3/2 * 5/4 * 7/6 = 2.1875).
Для потомков, вот переписывание прекрасного алгоритма Уилла Несса для Python 3. Необходимы некоторые изменения (итераторы больше не имеют .next() методы, но есть новый next() встроенная функция). Другие изменения для удовольствия (используя новый yield from <iterable> заменяет четыре yield заявления в оригинале. Больше для удобства чтения (я не фанат чрезмерного использования;-) 1-буквенные имена переменных).
Это значительно быстрее, чем оригинал, но не по алгоритмическим причинам. Ускорение в основном связано с удалением оригинала add() функция, делая это встроенным вместо этого.
def psieve():
import itertools
yield from (2, 3, 5, 7)
D = {}
ps = psieve()
next(ps)
p = next(ps)
assert p == 3
psq = p*p
for i in itertools.count(9, 2):
if i in D: # composite
step = D.pop(i)
elif i < psq: # prime
yield i
continue
else: # composite, = p*p
assert i == psq
step = 2*p
p = next(ps)
psq = p*p
i += step
while i in D:
i += step
D[i] = step
Изначально это не мой код, но стоит опубликовать. Оригинал можно найти здесь: http://code.activestate.com/recipes/117119/
def gen_primes():
D = {}
q = 2 # first integer to test for primality.
while True:
if q not in D:
# not marked composite, must be prime
yield q
#first multiple of q not already marked
D[q * q] = [q]
else:
for p in D[q]:
D.setdefault(p + q, []).append(p)
# no longer need D[q], free memory
del D[q]
q += 1
Это генератор, так что используйте его, как и любой другой.
primes = gen_primes()
for p in primes:
print p
Требуется 1,62 с, чтобы сгенерировать и поместить в набор 1 миллион простых чисел на моем рабочем столе.
Сделайте сегментированное сито, где размер сегмента определяется доступной памятью или максимальным размером набора битов.
Для каждого сегмента представляют числа в некотором интервале [n; n + сегмент_размер) в качестве набора битов и просеивания со всеми простыми числами ниже квадратного корня верхней границы.
Использование набора битов требует меньше памяти, чем хеш-таблица или древовидная структура данных, потому что вы работаете с плотными наборами чисел.
Это похоже, но примерно на 5% быстрее, чем erat2a.
import itertools as it
def erat2b( ):
D = { }
yield 2
for q in it.islice(it.count(3), 0, None, 2):
p = D.pop(q, None)
if p is None:
D[q*q] = q
yield q
else:
# Only this part is replaced.
q += p
while q in D:
q += p
D[q] = p
Еще один способ сделать это:
import itertools
def primeseq():
prime = [2]
num = 0
yield 2
for i in itertools.count(3, 2):
is_prime = True
for num in prime:
if i % num == 0:
is_prime = False
break
elif num ** 2 > i:
break
if is_prime:
prime.append(i)
yield i
И еще один ответ, более эффективный, чем мой erat3 ответ здесь:
import heapq
def heapprimegen():
hp= []
yield 2
yield 3
cn= 3
nn, inc= 3, 6
while 1:
while cn < nn:
yield cn
heapq.heappush(hp, (3*cn, 2*cn))
cn+= 2
cn= nn+2
nn, inc= heapq.heappushpop(hp, (nn+inc, inc))
Он поддерживает кучу (список) простых кратных, а не словарь. Очевидно, он теряет скорость.
Вот сложная реализация на основе кучи, которая не намного быстрее, чем другие реализации на основе кучи (см. Сравнение скорости в другом моем ответе), но использует намного меньше памяти.
Эта реализация использует две кучи (tu и wv), которые содержат одинаковые числовые элементы. Каждый элемент представляет собой пару int. Чтобы найти все простые числа до q**2 (где q простое число), каждая куча будет содержать не более 2*pi(q-1) элементы, где pi(x) число положительных простых чисел не больше x, Таким образом, общее число целых не более 4*pi(floor(sqrt(n))), (Мы могли бы увеличить коэффициент памяти на 2, поместив вдвое меньше содержимого в кучу, но это замедлило бы алгоритм.)
Другие подходы на основе dict и heap (например, erat2b, heap_prime_gen_squares и heapprimegen) выше хранят около 2*pi(n)'целых чисел, потому что они расширяют свою кучу или dict каждый раз, когда находят простое число. Для сравнения: чтобы найти простые числа 1_000_000, эта реализация хранит менее 4141 целых чисел, другие реализации хранят более 1_000_000 целых чисел.
import heapq
def heap_prime_gen_smallmem():
yield 2
yield 3
f = 5
fmar3 = 2
q = 7
q6 = 7 * 6
qmar3 = 4
tu = [(25, 30), (35, 30)]
vw = [(25, 30), (35, 30)]
while True:
qmar3 += 2
if qmar3 == 6:
qb = q + 4
q6b = q6 + 24
qmar3 = 2
else:
qb = q + 2
q6b = q6 + 12
if q < tu[0][0]:
d = q * q
while f < d:
a, b = vw[0]
if f < a:
yield f
else:
a, b = vw[0]
heapq.heapreplace(vw, (a + b, b))
a, b = vw[0]
while f >= a:
heapq.heapreplace(vw, (a + b, b))
a, b = vw[0]
fmar3 += 2
if fmar3 == 6:
f += 4
fmar3 = 2
else:
f += 2
c = q * qb
heapq.heappush(tu, (d, q6))
heapq.heappush(tu, (c, q6))
heapq.heappush(vw, (d, q6))
heapq.heappush(vw, (c, q6))
else:
a, b = tu[0]
heapq.heapreplace(tu, (a + b, b))
a, b = tu[0]
while q >= a:
heapq.heapreplace(tu, (a + b, b))
a, b = tu[0]
q = qb
q6 = q6b
Вот довольно быстрый бесконечный генератор, написанный на Python2, но легко адаптируемый к Python3. Чтобы использовать его для добавления простых чисел до 10**9, используйте следующее:
from itertools import takewhile
from functools import partial
from operator import gt
print (sum(takewhile(partial(gt, 10**9), prime_gen_inf())))
Это сегментированное сито, более быстрое, но явно менее элегантное, чем алгоритм Уилла Несса.
from operator import mul
from functools import reduce
def prod(x): return reduce(mul, x, 1)
def build_sieve(wheel):
w = prod(wheel)
w_phi = prod([p-1 for p in wheel])
rems = [a for a in range(w) if all(a % p for p in wheel)]
assert len(rems) == w_phi
inv = {a:pow(a, w_phi - 1, w) for a in rems}
try:
known_p = wheel + rems[1 : rems.index(rems[1]*rems[1])]
except ValueError:
known_p = wheel + rems[1:]
return wheel, w, w_phi, rems, inv, known_p
#Adjust the chunk variable based on your computer's architecture.
#
#Adjust the line with #! if you don't need "true" infinite. If you don't need
#primes larger than 1<<32, use array('H', []), if 1<<64 use 'L', if 1<<128 (in
#Python3) use 'Q', otherwise use empty list [].
#To save memory, comment out the lines with #*, and uncomment the commented-out
#lines
import itertools
from itertools import islice, count, compress, izip
chain_f = itertools.chain.from_iterable
from array import array
def prime_gen_inf(chunk=250000, sieve_info = build_sieve([2,3,5,7])):
""" Indefinitely yields primes """
wheel, w, w_phi, rems, inv, known_p = sieve_info
for p in known_p: yield p
new_n = 0;
while True:
size = min(chunk, (p * p - new_n) / w)
sieve = bytearray([1]) * size * w_phi
n, new_n = new_n, new_n + size * w
if not n:
zero = bytearray([0])
seen = len(known_p) - len(wheel) + 1
sieve[:seen:1] = zero * seen
p_gen = islice(prime_gen_inf(), len(wheel), None)
new_p = next(p_gen)
ps = [] #! array('H', [])
p_invs = bytearray([]) #*
while new_p * new_p < new_n:
ps.append(new_p)
p_invs.append(inv[new_p % w]) #*
new_p = next(p_gen)
for p, p_inv, modp in izip(ps, p_invs, [-n % p for p in ps]): #*
s = [(modp + p * (p_inv * (r - modp) % w)) / w for r in rems] #*
#for p in ps:
# s = [(-n%p + p * (inv[p%w] * (r - -n%p) % w)) / w for r in rems]
for i, start in enumerate(s):
slice_size = ((size - start - 1) / p + 1)
sieve[i + start * w_phi :: p * w_phi] = zero * slice_size
for p in compress(chain_f(izip(*[count(n+r, w) for r in rems])), sieve):
yield p
Вот генератор, который немного правдивее того, как это делается в Haskell: фильтрация по композициям известных простых чисел, затем добавление оставшихся простых чисел в список.
def gen_primes():
primes = []
i = 2
while True:
prime = True
for p in primes:
if not (i % p):
prime = False
break
if prime:
yield i
primes.append(i)
i += 1
Вот простой, но не очень медленный, использующий кучу вместо dict:
import heapq
def heap_prime_gen_squares():
yield 2
yield 3
h = [(9, 6)]
n = 5
while True:
a, b = h[0]
while n < a:
yield n
heapq.heappush(h, (n * n, n << 1))
n += 2
heapq.heapreplace(h, (a + b, b)) # Replace h[0], which is still (a, b).
Мои измерения скорости пользовательского времени на первые 1 миллион простых чисел (чем меньше число, тем лучше):
- postponed_sieve (на основе dict): 8,553
- erat2b (на основе dict): 9,513 с
- erat2a (на основе диктов): 10,313 с
- heap_prime_gen_smallmem (на основе кучи): 23,935 с
- heap_prime_gen_squares (на основе кучи): 27,302 с
- Heapprimegen (на основе dict): 145,029 с
Таким образом, основанные на диктате подходы кажутся самыми быстрыми.
Несколько раз назад я написал статью о бесконечном генераторе простых чисел:
http://stacktrace.it/2008/01/progetto-eulero-problema-3/
Это на итальянском, но вы можете сделать досадный перевод с помощью Google: http://tinyurl.com/yzpyeom
Я знаю, что этот пост старый, но я сам столкнулся с этим вопросом... Следующий код основан на очень простой идее: растущем сите Эратосфена. Это решение действительно медленнее, чем лучшие здесь, но его легко понять и оно предназначено для чтения...
Я использовал целые числа для хранения результатов сита. В двоичном формате целое число представляет собой список 0с и 1s, 0 в положении i если i не простое число, 1 если это может быть простое число. Необходимая бесконечность является результатом того, что целые числа Python 3 не ограничены.
def primes():
container, size = 1 << 2, 3 # we start with 0b100 (from right to left: 0 and 1 are not primes, 2 is
last_prime = 1
while True:
prime = next((j for j in range(last_prime+1, size) if container & 1 << j), None) # find the next prime
while not prime:
container, size = expand(container, size, 2**16) # add 65536 cells and sieve the container
prime = next((j for j in range(last_prime+1, size) if container & 1 << j), None)
yield prime
last_prime = prime
Как расширить контейнер? Просто добавьте кучу 1s слева от контейнера (в двоичном формате) и просеять их. Это идентично стандартному сита с небольшой разницей. В стандартном сите, если мы найдем простое число i, мы начинаем пересекать клетки в i*iс шагом i,
Здесь это может быть сделано для первой части контейнера. Нам просто нужно начать с начала новой части контейнера, если он дальше, чем i*i,
def expand(container, size, n):
new_size = size + n
container += (1 << (new_size + 1) - 1) - (1 << size) # add n 1's
for i in range(2, new_size):
if container & (1 << i): # i is a prime
t = sum(1 << j for j in range(max(i, size // i)*i, new_size, i)) # set 1 for all mutiple
container &= ~t # cross the cells
return container, new_size
Тест на миллион простых чисел:
import itertools
assert 78498 == len(list(itertools.takewhile(lambda p: p<1000000, primes())))
def generate_primes():
""" primes generator but not the most efficient solution"""
yield 2
n = 1
while True:
n += 2
if 0 not in (n % x for x in range(3, n//2, 2)):
yield n