Ответ 1

Мне все еще нравится то, что я написал здесь (рецепт поваренной книги со многими другими авторами) - он показывает, как нет сита Эратосфена внутренние ограничения, комментарии и обсуждение, я считаю, делают это совершенно ясно. Это было недавно обсуждено в разделе "Переполнение стека" (по-моему, поиск имен авторов), и кто-то предложил значительно более быструю (но IMHO менее ясную) версию; -).

Ответ 2

"Если я увижу дальше..."

Функция 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 элементов (имеется 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

Ответ 3

Поскольку 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).

Ответ 4

Для потомков здесь переписывается красивый алгоритм 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

Ответ 5

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

Ответ 6

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

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

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

Ответ 7

Это похоже, но на 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

Ответ 8

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

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

Ответ 9

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

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

Ответ 10

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

В этой реализации используются две кучи (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

Ответ 11

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

Ответ 12

Я написал статью о генераторе бесконечных простых чисел несколько раз назад:

http://stacktrace.it/2008/01/progetto-eulero-problema-3/

Это по-итальянски, но у вас может быть надоедливый перевод с помощью Google: http://tinyurl.com/yzpyeom

Ответ 13

Вот простой, но не очень медленный, используя кучу вместо 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.553s
• erat2b (на основе dict): 9.513s
• erat2a (на основе dict): 10.313s
• heap_prime_gen_smallmem (на основе кучи): 23.935s
• heap_prime_gen_squares (на основе кучи): 27.302s
• heapprimegen (на основе dict): 145.029s

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

Ответ 14

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