Как реализовать эффективный бесконечный генератор простых чисел в Python? - PullRequest
58 голосов
/ 06 февраля 2010

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

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

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

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

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

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

Спасибо, что читаете (и пишете;))!

Ответы [ 13 ]

69 голосов
/ 26 сентября 2010

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

Функция 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. Таким образом, после получения начальных трех простых чисел, мы начинаем с 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

Код теста

Модуль

A 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
65 голосов
/ 24 мая 2012

Поскольку 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) ). Следовательно, временная сложность также улучшается, то есть она работает быстрее .

Создает «скользящее сито» в качестве словаря текущих кратных значений каждого базового простого числа (т. Е. Ниже квадрата текущей производственной точки) вместе с их значениями 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 * Код на основе 1039 * работает ~ в 2,15 раза быстрее (что очень близко к теоретическому улучшению 3/2 * 5/4 * 7/6 = 2.1875).

37 голосов
/ 16 октября 2013

Для потомков, вот переписывание прекрасного алгоритма Уилла Несса для 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
5 голосов
/ 06 февраля 2010

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

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

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

5 голосов
/ 06 февраля 2010

Изначально это не мой код, но стоит опубликовать. Оригинал можно найти здесь: 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 миллион простых чисел на моем рабочем столе.

2 голосов
/ 01 августа 2012

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

Эта реализация использует две кучи (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
2 голосов
/ 05 февраля 2012

Еще один способ сделать это:

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
2 голосов
/ 05 декабря 2011

И еще один ответ, более эффективный по памяти, чем мой 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))

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

1 голос
/ 06 ноября 2015

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

Вот простой, но не очень медленный, использующий кучу вместо 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 (на основе диктов): 9,513 с
  • erat2a (на основе диктов): 10,313 с
  • heap_prime_gen_smallmem (на основе кучи): 23,935 с
  • heap_prime_gen_squares (на основе кучи): 27,302 с
  • Гепапримеген (на основе диктов): 145,029 с

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

...