«Если бы я видел дальше…»

Функцию 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% с помощью следующих изменений (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 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

Поскольку OP запрашивает эффективную реализацию, здесь существенное улучшение код активного состояния 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)). Следовательно, временная сложность также улучшается, т. Е. выполняется быстрее.

Создает "скользящее решето" как словарь текущих кратных каждому базовому простому числу (т.е. ниже sqrt текущей точки производства) вместе с их значениями step:

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

(старый исходный код здесь был отредактирован, чтобы включить изменения, указанные в ответе Тим Петерс, см. ниже). см. также this для обсуждения по теме.

Аналогичное колесо 2-3-5-7 - код на основе работает примерно в 2,15 раза быстрее (что очень близко к теоретическому улучшению 3/2 * 5/4 * 7/6 = 2.1875).

Для потомков - это переписанный вариант прекрасного алгоритма Уилла Несса для Python 3. Необходимы некоторые изменения (итераторы больше не имеют .next() методов , но есть новая встроенная функция next()). Другие изменения сделаны для развлечения (использование нового yield from <iterable> заменяет четыре оператора yield в оригинале. Больше для удобства чтения (я не фанат чрезмерного использования ;-) однобуквенных имен переменных).

Это значительно быстрее оригинала, но не по алгоритмам. Ускорение в основном связано с удалением исходной функции 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 миллиона простых чисел на моем рабочем столе, требуется 1,62 с.

Сделайте сегментированное сито, где размер сегмента определяется доступной памятью или максимальным размером битового набора.

Для каждого сегмента представляют числа в некотором интервале [n; n + segment_size) как набор бит и решето со всеми простыми числами ниже квадратного корня из верхней границы.

Использование набора битов требует меньше памяти, чем хэш-таблица или древовидная структура данных, потому что вы работаете с плотными наборами чисел.

Другой способ сделать это:

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))

Он поддерживает кучу (список) простых кратных, а не словарь. Очевидно, он немного теряет скорость.

Вот довольно быстрый бесконечный генератор, написанный на 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

Вот сложная реализация на основе кучи, которая не намного быстрее, чем другие реализации на основе кучи (см. Сравнение скорости в другом моем ответе), но она использует гораздо меньше памяти.

Эта реализация использует две кучи (tu и wv), которые содержат одинаковые числовые элементы. Каждый элемент представляет собой пару int. Чтобы найти все простые числа до q**2 (где q - простое число), каждая куча будет содержать не более 2*pi(q-1) элементов, где pi(x) - количество положительных простых чисел, не превышающих x. Таким образом, общее количество целых чисел не превышает 4*pi(floor(sqrt(n))). (Мы могли бы увеличить объем памяти вдвое, поместив в кучу вдвое меньше данных, но это замедлило бы алгоритм.)

Другие подходы на основе dict и кучи (например, 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

Вот простой, но не очень медленный вариант, использующий кучу вместо 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).

Мои измерения скорости пользовательского времени для первого миллиона простых чисел (меньшие числа лучше):

• postponed_sieve (на основе dict): 8.553 с
• erat2b (на основе диктовки): 9,513 с
• erat2a (на основе диктовки): 10,313 с
• heap_prime_gen_smallmem (на основе кучи): 23.935 с
• heap_prime_gen_squares (на основе кучи): 27,302 с
• heapprimegen (на основе диктовки): 145.029 с

Так что подходы, основанные на диктовке, кажутся самыми быстрыми.

Вот генератор, который немного ближе к тому, как это делается в 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

Я знаю, что сообщение старое, но я сам натолкнулся на этот вопрос ... Следующий код основан на очень простой идее: растущее сито Эратосфена. Это решение действительно медленнее, чем лучшие здесь, но его легко понять и оно разработано так, чтобы его можно было читать ...

Я использовал целые числа для хранения результатов сита. В двоичном формате целое число представляет собой список 0s и 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

Как расширить емкость? Просто добавьте кучу 1 слева от контейнера (в двоичном формате) и просейте их. Оно идентично стандартному сито с небольшой разницей. В стандартном сите, если мы находим простое число 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())))