Самый быстрый способ перечислить все простые числа ниже N - PullRequest
338 голосов
/ 15 января 2010

Это лучший алгоритм, который я мог придумать.

def get_primes(n):
    numbers = set(range(n, 1, -1))
    primes = []
    while numbers:
        p = numbers.pop()
        primes.append(p)
        numbers.difference_update(set(range(p*2, n+1, p)))
    return primes

>>> timeit.Timer(stmt='get_primes.get_primes(1000000)', setup='import   get_primes').timeit(1)
1.1499958793645562

Можно ли сделать это еще быстрее?

Этот код имеет недостаток: поскольку numbers является неупорядоченным набором, нет гарантии, что numbers.pop() удалит самый низкий номер из набора. Тем не менее, это работает (по крайней мере для меня) для некоторых входных чисел:

>>> sum(get_primes(2000000))
142913828922L
#That's the correct sum of all numbers below 2 million
>>> 529 in get_primes(1000)
False
>>> 529 in get_primes(530)
True

Ответы [ 31 ]

4 голосов
/ 21 января 2010

Если у вас есть контроль над N, самый быстрый способ перечислить все простые числа - это предварительно их вычислить. Шутки в сторону. Предварительные вычисления - это способ игнорирования оптимизации.

4 голосов
/ 21 января 2010

Детерминистическая реализация критерия Примитивности Миллера-Рабина в предположении, что N <9,080,191 </p>

import sys
import random

def miller_rabin_pass(a, n):
    d = n - 1
    s = 0
    while d % 2 == 0:
        d >>= 1
        s += 1

    a_to_power = pow(a, d, n)
    if a_to_power == 1:
        return True
    for i in xrange(s-1):
        if a_to_power == n - 1:
            return True
        a_to_power = (a_to_power * a_to_power) % n
    return a_to_power == n - 1


def miller_rabin(n):
    for a in [2, 3, 37, 73]:
      if not miller_rabin_pass(a, n):
        return False
    return True


n = int(sys.argv[1])
primes = [2]
for p in range(3,n,2):
  if miller_rabin(p):
    primes.append(p)
print len(primes)

Согласно статье в Википедии (http://en.wikipedia.org/wiki/Miller–Rabin_primality_test) тестирование N

И я адаптировал исходный код из вероятностной реализации оригинального теста Миллера-Рабина, найденного здесь: http://en.literateprograms.org/Miller-Rabin_primality_test_(Python)

4 голосов
/ 21 января 2010

Вот код, который я обычно использую для генерации простых чисел в Python:

$ python -mtimeit -s'import sieve' 'sieve.sieve(1000000)' 
10 loops, best of 3: 445 msec per loop
$ cat sieve.py
from math import sqrt

def sieve(size):
 prime=[True]*size
 rng=xrange
 limit=int(sqrt(size))

 for i in rng(3,limit+1,+2):
  if prime[i]:
   prime[i*i::+i]=[False]*len(prime[i*i::+i])

 return [2]+[i for i in rng(3,size,+2) if prime[i]]

if __name__=='__main__':
 print sieve(100)

Он не может конкурировать с более быстрыми решениями, размещенными здесь, но, по крайней мере, это чистый питон.

Спасибо за размещение этого вопроса. Я действительно многому научился сегодня.

3 голосов
/ 15 января 2010

Для самого быстрого кода лучшее решение - цифровое решение. Однако по чисто академическим причинам я публикую свою версию на чистом Python, которая чуть менее чем на 50% быстрее, чем версия поваренной книги, опубликованная выше. Поскольку я делаю весь список в памяти, вам нужно достаточно места, чтобы вместить все, но, похоже, он достаточно хорошо масштабируется.

def daniel_sieve_2(maxNumber):
    """
    Given a number, returns all numbers less than or equal to
    that number which are prime.
    """
    allNumbers = range(3, maxNumber+1, 2)
    for mIndex, number in enumerate(xrange(3, maxNumber+1, 2)):
        if allNumbers[mIndex] == 0:
            continue
        # now set all multiples to 0
        for index in xrange(mIndex+number, (maxNumber-3)/2+1, number):
            allNumbers[index] = 0
    return [2] + filter(lambda n: n!=0, allNumbers)

И результаты:

>>>mine = timeit.Timer("daniel_sieve_2(1000000)",
...                    "from sieves import daniel_sieve_2")
>>>prev = timeit.Timer("get_primes_erat(1000000)",
...                    "from sieves import get_primes_erat")
>>>print "Mine: {0:0.4f} ms".format(min(mine.repeat(3, 1))*1000)
Mine: 428.9446 ms
>>>print "Previous Best {0:0.4f} ms".format(min(prev.repeat(3, 1))*1000)
Previous Best 621.3581 ms
3 голосов
/ 05 октября 2016

Это все написано и проверено. Таким образом, нет необходимости изобретать велосипед.

python -m timeit -r10 -s"from sympy import sieve" "primes = list(sieve.primerange(1, 10**6))"

дает нам рекордный 12,2 мсек !

10 loops, best of 10: 12.2 msec per loop

Если это не достаточно быстро, вы можете попробовать PyPy:

pypy -m timeit -r10 -s"from sympy import sieve" "primes = list(sieve.primerange(1, 10**6))"

, что приводит к:

10 loops, best of 10: 2.03 msec per loop

Ответ с 247 списками голосов за 15,9 мс для лучшего решения. Сравните это !!!

3 голосов
/ 18 октября 2016

Я протестировал некоторые функции unutbu , я вычислил их с числом голодных миллионов

Победителями становятся функции, использующие библиотеку numpy,

Примечание : Было бы также интересно провести тест на использование памяти:)

Computation time result

Пример кода

Полный код в моем репозитории github

#!/usr/bin/env python

import lib
import timeit
import sys
import math
import datetime

import prettyplotlib as ppl
import numpy as np

import matplotlib.pyplot as plt
from prettyplotlib import brewer2mpl

primenumbers_gen = [
    'sieveOfEratosthenes',
    'ambi_sieve',
    'ambi_sieve_plain',
    'sundaram3',
    'sieve_wheel_30',
    'primesfrom3to',
    'primesfrom2to',
    'rwh_primes',
    'rwh_primes1',
    'rwh_primes2',
]

def human_format(num):
    # /551796/formatirovanie-dlinnyh-chisel-v-vide-strok-v-pythontab-top
    magnitude = 0
    while abs(num) >= 1000:
        magnitude += 1
        num /= 1000.0
    # add more suffixes if you need them
    return '%.2f%s' % (num, ['', 'K', 'M', 'G', 'T', 'P'][magnitude])


if __name__=='__main__':

    # Vars
    n = 10000000 # number itereration generator
    nbcol = 5 # For decompose prime number generator
    nb_benchloop = 3 # Eliminate false positive value during the test (bench average time)
    datetimeformat = '%Y-%m-%d %H:%M:%S.%f'
    config = 'from __main__ import n; import lib'
    primenumbers_gen = {
        'sieveOfEratosthenes': {'color': 'b'},
        'ambi_sieve': {'color': 'b'},
        'ambi_sieve_plain': {'color': 'b'},
         'sundaram3': {'color': 'b'},
        'sieve_wheel_30': {'color': 'b'},
# # #        'primesfrom2to': {'color': 'b'},
        'primesfrom3to': {'color': 'b'},
        # 'rwh_primes': {'color': 'b'},
        # 'rwh_primes1': {'color': 'b'},
        'rwh_primes2': {'color': 'b'},
    }


    # Get n in command line
    if len(sys.argv)>1:
        n = int(sys.argv[1])

    step = int(math.ceil(n / float(nbcol)))
    nbs = np.array([i * step for i in range(1, int(nbcol) + 1)])
    set2 = brewer2mpl.get_map('Paired', 'qualitative', 12).mpl_colors

    print datetime.datetime.now().strftime(datetimeformat)
    print("Compute prime number to %(n)s" % locals())
    print("")

    results = dict()
    for pgen in primenumbers_gen:
        results[pgen] = dict()
        benchtimes = list()
        for n in nbs:
            t = timeit.Timer("lib.%(pgen)s(n)" % locals(), setup=config)
            execute_times = t.repeat(repeat=nb_benchloop,number=1)
            benchtime = np.mean(execute_times)
            benchtimes.append(benchtime)
        results[pgen] = {'benchtimes':np.array(benchtimes)}

fig, ax = plt.subplots(1)
plt.ylabel('Computation time (in second)')
plt.xlabel('Numbers computed')
i = 0
for pgen in primenumbers_gen:

    bench = results[pgen]['benchtimes']
    avgs = np.divide(bench,nbs)
    avg = np.average(bench, weights=nbs)

    # Compute linear regression
    A = np.vstack([nbs, np.ones(len(nbs))]).T
    a, b = np.linalg.lstsq(A, nbs*avgs)[0]

    # Plot
    i += 1
    #label="%(pgen)s" % locals()
    #ppl.plot(nbs, nbs*avgs, label=label, lw=1, linestyle='--', color=set2[i % 12])
    label="%(pgen)s avg" % locals()
    ppl.plot(nbs, a * nbs + b, label=label, lw=2, color=set2[i % 12])
print datetime.datetime.now().strftime(datetimeformat)

ppl.legend(ax, loc='upper left', ncol=4)

# Change x axis label
ax.get_xaxis().get_major_formatter().set_scientific(False)
fig.canvas.draw()
labels = [human_format(int(item.get_text())) for item in ax.get_xticklabels()]

ax.set_xticklabels(labels)
ax = plt.gca()

plt.show()
3 голосов
/ 03 сентября 2010

Немного другая реализация полусита с использованием Numpy:

http://rebrained.com/?p=458

import math
import numpy
def prime6(upto):
    primes=numpy.arange(3,upto+1,2)
    isprime=numpy.ones((upto-1)/2,dtype=bool)
    for factor in primes[:int(math.sqrt(upto))]:
        if isprime[(factor-2)/2]: isprime[(factor*3-2)/2:(upto-1)/2:factor]=0
    return numpy.insert(primes[isprime],0,2)

Может кто-нибудь сравнить это с другими временами? На моей машине это выглядит довольно сравнимо с другим половинным ситом Numpy.

3 голосов
/ 19 июля 2017

для Python 3

def rwh_primes2(n):
    correction = (n%6>1)
    n = {0:n,1:n-1,2:n+4,3:n+3,4:n+2,5:n+1}[n%6]
    sieve = [True] * (n//3)
    sieve[0] = False
    for i in range(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+4*k-2*k*(i&1))//3::2*k]=[False]*((n//6-(k*k+4*k-2*k*(i&1))//6-1)//k+1)
    return [2,3] + [3*i+1|1 for i in range(1,n//3-correction) if sieve[i]]
2 голосов
/ 22 сентября 2013

Я знаю, что конкурс закрыт на несколько лет. ...

Тем не менее, это мое предложение для простого сита Python, основанное на пропуске кратных 2, 3 и 5 с использованием соответствующих шагов при обработке сита вперед. Тем не менее, для N <10 ^ 9 он на самом деле медленнее, чем у @Robert William Hanks превосходных решений rwh_primes2 и rwh_primes1. Использование ситового массива ctypes.c_ushort выше 1,5 * 10 ^ 8 позволяет адаптировать его к ограничениям памяти. </p>

10 ^ 6

$ python -mtimeit -s "import primeSieveSpeedComp" "primeSieveSpeedComp.primeSieveSeq (1000000)" 10 циклов, лучшее из 3: 46,7 мсек на цикл

для сравнения: $ python -mtimeit -s "import primeSieveSpeedComp" "primeSieveSpeedComp.rwh_primes1 (1000000)" 10 циклов, лучшее из 3: 43.2 мсек на цикл для сравнения: $ python -m timeit -s "import primeSieveSpeedComp" «primeSieveSpeedComp.rwh_primes2 (1000000)» 10 циклов, лучшее из 3: 34,5 мсек на цикл

10 ^ 7

$ python -mtimeit -s "import primeSieveSpeedComp" "primeSieveSpeedComp.primeSieveSeq (10000000)" 10 циклов, лучше всего 3: 530 мсек на цикл

для сравнения: $ python -mtimeit -s "import primeSieveSpeedComp" "primeSieveSpeedComp.rwh_primes1 (10000000)" 10 циклов, лучшее из 3: 494 мсек на цикл для сравнения: $ python -m timeit -s "import primeSieveSpeedComp" "primeSieveSpeedComp.rwh_primes2 (10000000)" 10 циклов, лучшее из 3: 375 мсек на цикл

10 ^ 8

$ python -mtimeit -s "import primeSieveSpeedComp" "primeSieveSpeedComp.primeSieveSeq (100000000)" 10 циклов, лучшее из 3: 5,55 с на цикл

для сравнения: $ python -mtimeit -s "import primeSieveSpeedComp" «primeSieveSpeedComp.rwh_primes1 (100000000)» 10 циклов, лучшее из 3: 5,33 сек за цикл для сравнения: $ python -m timeit -s "import primeSieveSpeedComp" "primeSieveSpeedComp.rwh_primes2 (100000000)" 10 циклов, лучшее из 3: 3,95 сек за цикл

10 ^ 9

$ python -mtimeit -s "import primeSieveSpeedComp" "primeSieveSpeedComp.primeSieveSeq (1000000000)" 10 циклов, лучшее из 3: 61,2 сек за цикл

для сравнения: $ python -mtimeit -n 3 -s "import primeSieveSpeedComp" "primeSieveSpeedComp.rwh_primes1 (1000000000)" 3 цикла, лучшее из 3: 97.8 сек за цикл

для сравнения: $ python -m timeit -s "import primeSieveSpeedComp" "primeSieveSpeedComp.rwh_primes2 (1000000000)" 10 циклов, лучшее из 3: 41,9 с за цикл

Вы можете скопировать приведенный ниже код в ubuntus primeSieveSpeedComp для просмотра этих тестов.

def primeSieveSeq(MAX_Int):
    if MAX_Int > 5*10**8:
        import ctypes
        int16Array = ctypes.c_ushort * (MAX_Int >> 1)
        sieve = int16Array()
        #print 'uses ctypes "unsigned short int Array"'
    else:
        sieve = (MAX_Int >> 1) * [False]
        #print 'uses python list() of long long int'
    if MAX_Int < 10**8:
        sieve[4::3] = [True]*((MAX_Int - 8)/6+1)
        sieve[12::5] = [True]*((MAX_Int - 24)/10+1)
    r = [2, 3, 5]
    n = 0
    for i in xrange(int(MAX_Int**0.5)/30+1):
        n += 3
        if not sieve[n]:
            n2 = (n << 1) + 1
            r.append(n2)
            n2q = (n2**2) >> 1
            sieve[n2q::n2] = [True]*(((MAX_Int >> 1) - n2q - 1) / n2 + 1)
        n += 2
        if not sieve[n]:
            n2 = (n << 1) + 1
            r.append(n2)
            n2q = (n2**2) >> 1
            sieve[n2q::n2] = [True]*(((MAX_Int >> 1) - n2q - 1) / n2 + 1)
        n += 1
        if not sieve[n]:
            n2 = (n << 1) + 1
            r.append(n2)
            n2q = (n2**2) >> 1
            sieve[n2q::n2] = [True]*(((MAX_Int >> 1) - n2q - 1) / n2 + 1)
        n += 2
        if not sieve[n]:
            n2 = (n << 1) + 1
            r.append(n2)
            n2q = (n2**2) >> 1
            sieve[n2q::n2] = [True]*(((MAX_Int >> 1) - n2q - 1) / n2 + 1)
        n += 1
        if not sieve[n]:
            n2 = (n << 1) + 1
            r.append(n2)
            n2q = (n2**2) >> 1
            sieve[n2q::n2] = [True]*(((MAX_Int >> 1) - n2q - 1) / n2 + 1)
        n += 2
        if not sieve[n]:
            n2 = (n << 1) + 1
            r.append(n2)
            n2q = (n2**2) >> 1
            sieve[n2q::n2] = [True]*(((MAX_Int >> 1) - n2q - 1) / n2 + 1)
        n += 3
        if not sieve[n]:
            n2 = (n << 1) + 1
            r.append(n2)
            n2q = (n2**2) >> 1
            sieve[n2q::n2] = [True]*(((MAX_Int >> 1) - n2q - 1) / n2 + 1)
        n += 1
        if not sieve[n]:
            n2 = (n << 1) + 1
            r.append(n2)
            n2q = (n2**2) >> 1
            sieve[n2q::n2] = [True]*(((MAX_Int >> 1) - n2q - 1) / n2 + 1)
    if MAX_Int < 10**8:
        return [2, 3, 5]+[(p << 1) + 1 for p in [n for n in xrange(3, MAX_Int >> 1) if not sieve[n]]]
    n = n >> 1
    try:
        for i in xrange((MAX_Int-2*n)/30 + 1):
            n += 3
            if not sieve[n]:
                r.append((n << 1) + 1)
            n += 2
            if not sieve[n]:
                r.append((n << 1) + 1)
            n += 1
            if not sieve[n]:
                r.append((n << 1) + 1)
            n += 2
            if not sieve[n]:
                r.append((n << 1) + 1)
            n += 1
            if not sieve[n]:
                r.append((n << 1) + 1)
            n += 2
            if not sieve[n]:
                r.append((n << 1) + 1)
            n += 3
            if not sieve[n]:
                r.append((n << 1) + 1)
            n += 1
            if not sieve[n]:
                r.append((n << 1) + 1)
    except:
        pass
    return r
2 голосов
/ 20 февраля 2013

Впервые использую python, поэтому некоторые методы, которые я использую в этом, могут показаться немного громоздкими. Я просто преобразовал свой код C ++ в Python, и это то, что у меня есть (хотя и немного медленный в Python)

#!/usr/bin/env python
import time

def GetPrimes(n):

    Sieve = [1 for x in xrange(n)]

    Done = False
    w = 3

    while not Done:

        for q in xrange (3, n, 2):
            Prod = w*q
            if Prod < n:
                Sieve[Prod] = 0
            else:
                break

        if w > (n/2):
            Done = True
        w += 2

    return Sieve



start = time.clock()

d = 10000000
Primes = GetPrimes(d)

count = 1 #This is for 2

for x in xrange (3, d, 2):
    if Primes[x]:
        count+=1

elapsed = (time.clock() - start)
print "\nFound", count, "primes in", elapsed, "seconds!\n"

pythonw Primes.py

Найдено 664579 простых чисел за 12,799119 секунд!

#!/usr/bin/env python
import time

def GetPrimes2(n):

    Sieve = [1 for x in xrange(n)]

    for q in xrange (3, n, 2):
        k = q
        for y in xrange(k*3, n, k*2):
            Sieve[y] = 0

    return Sieve



start = time.clock()

d = 10000000
Primes = GetPrimes2(d)

count = 1 #This is for 2

for x in xrange (3, d, 2):
    if Primes[x]:
        count+=1

elapsed = (time.clock() - start)
print "\nFound", count, "primes in", elapsed, "seconds!\n"

pythonw Primes2.py

Найдено 664579 простых чисел за 10.230172 секунды!

#!/usr/bin/env python
import time

def GetPrimes3(n):

    Sieve = [1 for x in xrange(n)]

    for q in xrange (3, n, 2):
        k = q
        for y in xrange(k*k, n, k << 1):
            Sieve[y] = 0

    return Sieve



start = time.clock()

d = 10000000
Primes = GetPrimes3(d)

count = 1 #This is for 2

for x in xrange (3, d, 2):
    if Primes[x]:
        count+=1

elapsed = (time.clock() - start)
print "\nFound", count, "primes in", elapsed, "seconds!\n"

python Primes2.py

Найдено 664579 простых чисел за 7,113776 секунд!

...