Это не домашнее задание, мне просто интересно.
INFINITE - это ключевое слово здесь.
Я хочу использовать его как для p в простых(). Я считаю, что это встроенная функция в Haskell.
Итак, ответ не может быть таким наивным, как "Просто сделать сито".
Прежде всего, вы не знаете, сколько последовательных простых чисел будет потреблено. Ну, предположим, что вы могли собрать 100 из них за раз. Используете ли вы такой же подход Сите, а также частоту формулы простых чисел?
Я предпочитаю неконкурентный подход.
Спасибо за чтение (и письмо;))!
Мне все еще нравится то, что я написал здесь (рецепт поваренной книги со многими другими авторами) - он показывает, как нет сита Эратосфена внутренние ограничения, комментарии и обсуждение, я считаю, делают это совершенно ясно. Это было недавно обсуждено в разделе "Переполнение стека" (по-моему, поиск имен авторов), и кто-то предложил значительно более быструю (но IMHO менее ясную) версию; -).
Функция erat2 из поваренной книги может быть дополнительно ускорена (примерно на 20-25%):
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, половина шагов избегается вместе с тестом для странности.
Если вы не возражаете против небольшого дополнительного причудливости, erat2 можно увеличить на 35-40% со следующими изменениями (NB: требуется 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. Таким образом, после получения начальных трех простых чисел, мы начинаем с 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() могут быть только такими -false-кандидатами.
На сервере 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
A primegen.py содержит функции erat2, erat2a и erat3. Здесь следует тестирование script:
#!/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
Поскольку OP запрашивает эффективную реализацию, здесь существенное улучшение в коде 2002 активного состояния от Дэвида Эппштейна /Alex Martelli (см. здесь в его ответ): не записывать основную информацию в словаре, пока ее квадрат не будет замечен среди кандидатов. Сложность пространства ниже до O (sqrt (n)) вместо O (n), для генерируемых n чисел ( π (sqrt (n log n)) ~ 2 sqrt (n log n)/log (n log n) ~ 2 sqrt (n/log n)). Следовательно, временная сложность также улучшается, т.е. работает быстрее.
Создает "скользящее сито" в качестве словаря текущих кратных значений каждого базового простого (т.е. ниже sqrt текущей производственной точки) вместе со своими значениями шага:
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 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 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.15x быстрее (что очень близко к теоретическому улучшению 3/2 * 5/4 * 7/6 = 2.1875).
Для потомков здесь переписывается красивый алгоритм Will Ness для Python 3. Требуются некоторые изменения (итераторы больше не имеют методов .next(), но там есть новая встроенная функция next()). Другие изменения для удовольствия (использование нового yield from <iterable> заменяет четыре оператора yield в оригинале. Больше для читаемости (я не поклонник чрезмерного использования;-) 1-буквенных имен переменных).
Это значительно быстрее, чем оригинал, но не по алгоритмическим причинам. Ускорение в основном связано с удалением исходной функции add(), делая это inline вместо.
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 миллион простых чисел, на моем рабочем столе требуется 1,62.
Сделайте сегментированное сито, где размер сегмента определяется доступной памятью или максимальным размером битового набора.
Для каждого сегмента представлены числа в некотором интервале [n; n + segment_size) в виде набора бит и сита со всеми простыми числами ниже квадратного корня верхней границы.
Использование набора бит использует меньше памяти, чем хеш-таблица или древовидная структура данных, потому что вы работаете с плотными наборами чисел.
Это похоже, но на 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
И еще один ответ, более эффективный с точки зрения памяти, чем мой 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))
Он поддерживает кучу (список) простых кратных, а не словарь. Очевидно, он теряет скорость.
Другой способ сделать это:
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
Вот сложная реализация на основе кучи, которая не намного быстрее, чем другие реализации на основе кучи (см. сравнение скорости в другом ответе), но она использует гораздо меньше памяти.
В этой реализации используются две кучи (tu и wv), содержащие одни и те же элементы числа. Каждый элемент представляет собой пару int. Чтобы найти все простые числа до q**2 (где q является простым), каждая куча будет содержать не более 2*pi(q-1) элементов, где pi(x) - число положительных простых чисел, не больших x. Таким образом, общее число целых чисел не более 4*pi(floor(sqrt(n))). (Мы могли бы получить коэффициент на 2 в памяти, нажимая половину всего материала в кучу, но это сделает алгоритм более медленным.)
Другие подходы, основанные на методе и куче (например, 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
Здесь генератор, немного верный тому, как это делается в 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
Я написал статью о генераторе бесконечных простых чисел несколько раз назад:
http://stacktrace.it/2008/01/progetto-eulero-problema-3/
Это по-итальянски, но у вас может быть надоедливый перевод с помощью Google: http://tinyurl.com/yzpyeom
Вот простой, но не очень медленный, используя кучу вместо 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 миллиона простых чисел (меньшие числа лучше):
Таким образом, подходы, основанные на dict, выглядят как самые быстрые.
Здесь довольно быстрый бесконечный генератор, написанный на 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 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