Быстрые способы нахождения всех простых чисел меньше N

#python #python_3x #алгоритм #оптимизация

Какие существуют алгоритмы наиболее быстрого нахождения всех простых чисел меньше
N кроме всем известных:


Решето Эратосфена
Решето Аткина
Решето Сундарама


И какие существуют быстрые реализации на Python (Vanilla, Numpy, etc.)?

Наивная реализация решета Эратосфена:

def eratosthenes2(n):
    multiples = set()
    for i in range(2, n+1):
        if i not in multiples:
            yield i
            multiples.update(range(i*i, n+1, i))


Аналогичный вопрос в англоязычной версии SO

PS идея создания данного вопроса связана с постоянно появляющимися новыми вопросами
об эффективной реализации нахождения всех простых чисел меньше N.
    

Ответы

Ответ 1

В качестве ответа я решил сравнить самые удачные реализации из ответов на данный
вопрос и сделать замеры времени (для этого я использовал модуль reporttime из этого
ответа @jfs:

primes.py



import numpy
from math import sqrt, ceil


def rwh_primes(n):
    """ Returns  a list of primes < n """
    # (c) Robert William Hanks - https://stackoverflow.com/a/3035188/5741205
    sieve = [True] * n
    for i in range(3,int(n**0.5)+1,2):
        if sieve[i]:
            sieve[i*i::2*i]=[False]*((n-i*i-1)//(2*i)+1)
    return [2] + [i for i in range(3,n,2) if sieve[i]]


def rwh_primes1(n):
    """ Returns  a list of primes < n """
    # (c) Robert William Hanks - https://stackoverflow.com/a/3035188/5741205
    sieve = [True] * (n//2)
    for i in range(3,int(n**0.5)+1,2):
        if sieve[i//2]:
            sieve[i*i//2::i] = [False] * ((n-i*i-1)//(2*i)+1)
    return [2] + [2*i+1 for i in range(1,n//2) if sieve[i]]


def rwh_primes2(n):
    """ Input n>=6, Returns a list of primes, 2 <= p < n """
    # (c) Robert William Hanks - https://stackoverflow.com/a/3035188/5741205
    n, correction = n-n%6+6, 2-(n%6>1)
    sieve = [True] * (n//3)
    for i in range(1,int(n**0.5)//3+1):
      if sieve[i]:
        k=3*i+1|1
        sieve[      k*k//3      ::2*k] = [False] * ((n//6-k*k//6-1)//k+1)
        sieve[k*(k-2*(i&1)+4)//3::2*k] = [False] * ((n//6-k*(k-2*(i&1)+4)//6-1)//k+1)
    return [2,3] + [3*i+1|1 for i in range(1,n//3-correction) if sieve[i]]


def primesfrom3to(n):
    """ Returns a array of primes, 3 <= p < n """
    # (c) Robert William Hanks - https://stackoverflow.com/a/3035188/5741205
    sieve = numpy.ones(n//2, dtype=numpy.bool)
    for i in range(3,int(n**0.5)+1,2):
        if sieve[i//2]:
            sieve[i*i//2::i] = False
    return 2*numpy.nonzero(sieve)[0][1::]+1


def primesfrom2to(n):
    """ Input n>=6, Returns a array of primes, 2 <= p < n """
    # (c) Robert William Hanks - https://stackoverflow.com/a/3035188/5741205
    sieve = numpy.ones(n//3 + (n%6==2), dtype=numpy.bool)
    for i in range(1,int(n**0.5)//3+1):
        if sieve[i]:
            k=3*i+1|1
            sieve[       k*k//3     ::2*k] = False
            sieve[k*(k-2*(i&1)+4)//3::2*k] = False
    return numpy.r_[2,3,((3*numpy.nonzero(sieve)[0][1:]+1)|1)]


def sieveOfEratosthenes(n):
    """sieveOfEratosthenes(n): return the list of the primes < n."""
    # Code from: , Nov 30 2006
    # http://groups.google.com/group/comp.lang.python/msg/f1f10ced88c68c2d
    if n <= 2:
        return []
    sieve = list(range(3, n, 2))
    top = len(sieve)
    for si in sieve:
        if si:
            bottom = (si*si - 3) // 2
            if bottom >= top:
                break
            sieve[bottom::si] = [0] * -((bottom - top) // si)
    return [2] + [el for el in sieve if el]


def sieveOfAtkin(end):
    """sieveOfAtkin(end): return a list of all the prime numbers , improved
    # Code: https://web.archive.org/web/20080324064651/http://krenzel.info/?p=83
    # Info: http://en.wikipedia.org/wiki/Sieve_of_Atkin
    assert end > 0
    lng = ((end-1) // 2)
    sieve = [False] * (lng + 1)

    x_max, x2, xd = int(sqrt((end-1)/4.0)), 0, 4
    for xd in range(4, 8*x_max + 2, 8):
        x2 += xd
        y_max = int(sqrt(end-x2))
        n, n_diff = x2 + y_max*y_max, (y_max << 1) - 1
        if not (n & 1):
            n -= n_diff
            n_diff -= 2
        for d in range((n_diff - 1) << 1, -1, -8):
            m = n % 12
            if m == 1 or m == 5:
                m = n >> 1
                sieve[m] = not sieve[m]
            n -= d

    x_max, x2, xd = int(sqrt((end-1) / 3.0)), 0, 3
    for xd in range(3, 6 * x_max + 2, 6):
        x2 += xd
        y_max = int(sqrt(end-x2))
        n, n_diff = x2 + y_max*y_max, (y_max << 1) - 1
        if not(n & 1):
            n -= n_diff
            n_diff -= 2
        for d in range((n_diff - 1) << 1, -1, -8):
            if n % 12 == 7:
                m = n >> 1
                sieve[m] = not sieve[m]
            n -= d

    x_max, y_min, x2, xd = int((2 + sqrt(4-8*(1-end)))/4), -1, 0, 3
    for x in range(1, x_max + 1):
        x2 += xd
        xd += 6
        if x2 >= end: y_min = (((int(ceil(sqrt(x2 - end))) - 1) << 1) - 2) << 1
        n, n_diff = ((x*x + x) << 1) - 1, (((x-1) << 1) - 2) << 1
        for d in range(n_diff, y_min, -8):
            if n % 12 == 11:
                m = n >> 1
                sieve[m] = not sieve[m]
            n += d

    primes = [2, 3]
    if end <= 3:
        return primes[:max(0,end-2)]

    for n in range(5 >> 1, (int(sqrt(end))+1) >> 1):
        if sieve[n]:
            primes.append((n << 1) + 1)
            aux = (n << 1) + 1
            aux *= aux
            for k in range(aux, end, 2 * aux):
                sieve[k >> 1] = False

    s  = int(sqrt(end)) + 1
    if s  % 2 == 0:
        s += 1
    primes.extend([i for i in range(s, end, 2) if sieve[i >> 1]])

    return primes


def sundaram3(max_n):
    # https://stackoverflow.com/questions/2068372/fastest-way-to-list-all-primes-below-n-in-python/2073279#2073279
    numbers = list(range(3, max_n+1, 2))
    half = (max_n)//2
    initial = 4

    for step in range(3, max_n+1, 2):
        for i in range(initial, half, step):
            numbers[i-1] = 0
        initial += 2*(step+1)

        if initial > half:
            return [2] + list(filter(None, numbers))




Замеры времени:

import os
import pandas as pd

os.chdir(r'C:\work\project\primes_timing')

import primes
from reporttime import get_functions_with_prefix, measure

funcs = get_functions_with_prefix('', module=primes)

d = {}
for n in [10**3, 10**4, 10**5, 10**6]:
    m = measure(funcs, args=[n])
    d[n] = pd.Series({b:a for a,b in m}) 
r = pd.DataFrame(d)


Вывод на печать:

name                     time ratio comment
primesfrom3to         21 usec  1.00 [1000]
sieveOfEratosthenes 40.1 usec  1.91 [1000]
rwh_primes          50.2 usec  2.39 [1000]
primesfrom2to       63.6 usec  3.02 [1000]
sundaram3           91.8 usec  4.37 [1000]
rwh_primes1          110 usec  5.25 [1000]
rwh_primes2          129 usec  6.12 [1000]
sieveOfAtkin         286 usec 13.61 [1000]

name                     time ratio comment
primesfrom3to       69.3 usec  1.00 [10000]
primesfrom2to        172 usec  2.48 [10000]
sieveOfEratosthenes  375 usec  5.41 [10000]
rwh_primes           447 usec  6.45 [10000]
rwh_primes2          503 usec  7.26 [10000]
rwh_primes1          577 usec  8.33 [10000]
sundaram3           1.22 msec 17.66 [10000]
sieveOfAtkin         2.2 msec 31.72 [10000]

name                     time ratio comment
primesfrom3to        306 usec  1.00 [100000]
primesfrom2to        371 usec  1.21 [100000]
sieveOfEratosthenes 3.82 msec 12.47 [100000]
rwh_primes          3.94 msec 12.84 [100000]
rwh_primes2         4.19 msec 13.68 [100000]
rwh_primes1         5.37 msec 17.53 [100000]
sundaram3           15.2 msec 49.57 [100000]
sieveOfAtkin        21.3 msec 69.47 [100000]

name                     time  ratio comment
primesfrom2to       2.29 msec   1.00 [1000000]
primesfrom3to       3.05 msec   1.33 [1000000]
rwh_primes2         38.7 msec  16.91 [1000000]
rwh_primes1         48.3 msec  21.12 [1000000]
rwh_primes          52.5 msec  22.96 [1000000]
sieveOfEratosthenes   61 msec  26.64 [1000000]
sieveOfAtkin         205 msec  89.64 [1000000]
sundaram3            245 msec 107.27 [1000000]


DataFrame с результатами:

In [67]: r
Out[67]:
                      1000      10000     100000    1000000
primesfrom2to        0.000064  0.000172  0.000371  0.002288
primesfrom3to        0.000021  0.000069  0.000306  0.003046
rwh_primes           0.000050  0.000447  0.003936  0.052516
rwh_primes1          0.000110  0.000577  0.005373  0.048307
rwh_primes2          0.000129  0.000503  0.004191  0.038678
sieveOfAtkin         0.000286  0.002197  0.021293  0.205074
sieveOfEratosthenes  0.000040  0.000375  0.003823  0.060955
sundaram3            0.000092  0.001223  0.015194  0.245407


график:

import matplotlib.pyplot as plt
import matplotlib

matplotlib.style.use('ggplot')

r.T.plot(loglog=True, grid=True, figsize=(8, 6))
plt.xlabel('N - generate all primes below N')
plt.ylabel('Time in seconds')
plt.show()


Ответ 2

Разбор работы самой быстрой из известных мне реализаций решета Эратосфена на Vanilla
(без дополнительных модулей) Python (c) Bruno Astrolino E Silva - я лишь добавил отладочную
информацию:

def primes(n):
    """ Returns  a list of primes < n """
    # (c) Robert William Hanks - https://stackoverflow.com/a/3035188/5741205
    sieve = [True] * n
    print("все чётные числа игнорируются и будут пропущены при возврате...\n")
    for i in range(3,int(n**0.5)+1,2):
        if sieve[i]:
            print('содержимое решета:\t{}'.format([x for x in range(3,n,2) if sieve[x]]))
            print(f'i:{i} вычёркиваем все числа кратные "{i}", начиная с "{i}^2":
{list(range(i*i, n, 2*i))}')
            sieve[i*i::2*i]=[False]*((n-i*i-1)//(2*i)+1)
            print(f'sieve[{i}*{i}::2*{i}]=[False]*(({n-i}*{i-1})//(2*{i})+1)')
            print('содержимое решета:\t{}'.format([x for x in range(3,n,2) if sieve[x]]))
            print('*' * 60)
    return [2] + [i for i in range(3,n,2) if sieve[i]]


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

Вывод для n=50:

In [165]: primes(50)
все чётные числа игнорируются и будут пропущены при возврате...

содержимое решета:      [3, 5, 7, 9, 11, 13, 15, 17, 19, 21, 23, 25, 27, 29, 31,
33, 35, 37, 39, 41, 43, 45, 47, 49]
i:3 вычёркиваем все числа кратные "3", начиная с "3^2": [9, 15, 21, 27, 33, 39, 45]
sieve[3*3::2*3]=[False]*((47*2)//(2*3)+1)
содержимое решета:      [3, 5, 7, 11, 13, 17, 19, 23, 25, 29, 31, 35, 37, 41, 43, 47, 49]
************************************************************
содержимое решета:      [3, 5, 7, 11, 13, 17, 19, 23, 25, 29, 31, 35, 37, 41, 43, 47, 49]
i:5 вычёркиваем все числа кратные "5", начиная с "5^2": [25, 35, 45]
sieve[5*5::2*5]=[False]*((45*4)//(2*5)+1)
содержимое решета:      [3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 49]
************************************************************
содержимое решета:      [3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 49]
i:7 вычёркиваем все числа кратные "7", начиная с "7^2": [49]
sieve[7*7::2*7]=[False]*((43*6)//(2*7)+1)
содержимое решета:      [3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47]
************************************************************
Out[165]: [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47