How to implement an efficient infinite generator of prime numbers in Python?

Это не домашнее задание, мне просто любопытно.

Ключевое слово здесь - БЕСКОНЕЧНОСТЬ.

Я хочу использовать это как for p in primes(). Я считаю, что это встроенная функция в Haskell.

Итак, ответ не может быть таким наивным, как "Просто сделай сито".

Прежде всего, вы не знаете, сколько последовательных простых чисел будет потреблено. Ну, предположим, вы могли бы придумать 100 из них за раз. Вы бы использовали тот же подход сита, а также формулу частоты простых чисел?

Я предпочитаю непараллельный подход.

Спасибо за чтение (и написание ;))!

“Если бы я видел дальше ...”

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 элементов (в каждых 30 числах 15 нечетных чисел, как выбрано 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 запрашивает эффективную реализацию, вот значительное улучшение кода active state 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 текущей производственной точки) вместе с их значениями шага:

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

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

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

Обновление 2022 года: Недавно я случайно наткнулся на эту старую штуку "NESL" 1990-х годов, которая на самом деле sqrt трюк с рекурсией, что и на Python.. Итак, ничто не ново под солнцем. :)

Код может быть напрямую дополнен, чтобы начать перечисление простых чисел непосредственно с заданного значения. Это можно увидеть в этой записи на основе JS.

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