Самый быстрый способ перечислить все простые числа ниже N

358

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

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
jbochi
источник
Фрагмент кода, о котором идет речь, гораздо быстрее, если числа, объявленные как числа = набор (диапазон (n, 2, -2)). Но не может победить sundaram3. Спасибо за вопрос.
Шехар
3
Было бы неплохо, если бы в ответах были версии функций Python 3.
Майкл Фукаракис
Конечно , есть библиотека , чтобы сделать это так , мы не должны свернуть наш собственный> XKCD обещал Python так просто import antigravity. Разве нет ничего подобного require 'prime'; Prime.take(10)(Руби)?
полковник Паник
2
@ColonelPanic Как это случилось, я обновил github.com/jaredks/pyprimesieve для Py3 и добавил в PyPi. Это, конечно, быстрее, чем эти, но не на несколько порядков - больше похоже на ~ 5 раз быстрее, чем лучшие версии.
Джаред
3
@ColonelPanic: Я думаю, что редактировать старые ответы, чтобы отметить, что они состарены, уместно, так как это делает его более полезным ресурсом. Если «принятый» ответ больше не является лучшим, возможно, отредактируйте примечание к вопросу с обновлением 2015 года, чтобы указать людям на текущий лучший метод.
Питер Кордес

Ответы:

366

Предупреждение: timeit результаты могут отличаться из-за различий в оборудовании или версии Python.

Ниже приведен скрипт, который сравнивает ряд реализаций:

Большое спасибо Стефану за то, что он привлек мое внимание к sieve_wheel_30. Заслуга Роберт Уильям Хэнкса для primesfrom2to, primesfrom3to, rwh_primes, rwh_primes1 и rwh_primes2.

Из простых методов Python, протестированных с psyco , для n = 1000000, rwh_primes1 был самым быстрым из протестированных.

+---------------------+-------+
| Method              | ms    |
+---------------------+-------+
| rwh_primes1         | 43.0  |
| sieveOfAtkin        | 46.4  |
| rwh_primes          | 57.4  |
| sieve_wheel_30      | 63.0  |
| rwh_primes2         | 67.8  |    
| sieveOfEratosthenes | 147.0 |
| ambi_sieve_plain    | 152.0 |
| sundaram3           | 194.0 |
+---------------------+-------+

Из простых методов Python, протестированных без psyco , для n = 1000000, rwh_primes2 был самым быстрым.

+---------------------+-------+
| Method              | ms    |
+---------------------+-------+
| rwh_primes2         | 68.1  |
| rwh_primes1         | 93.7  |
| rwh_primes          | 94.6  |
| sieve_wheel_30      | 97.4  |
| sieveOfEratosthenes | 178.0 |
| ambi_sieve_plain    | 286.0 |
| sieveOfAtkin        | 314.0 |
| sundaram3           | 416.0 |
+---------------------+-------+

Из всех протестированных методов, допускающих numpy , для n = 1000000, primesfrom2to был самым быстрым из протестированных.

+---------------------+-------+
| Method              | ms    |
+---------------------+-------+
| primesfrom2to       | 15.9  |
| primesfrom3to       | 18.4  |
| ambi_sieve          | 29.3  |
+---------------------+-------+

Время измерялось с помощью команды:

python -mtimeit -s"import primes" "primes.{method}(1000000)"

с {method}заменой каждого из имен методов.

primes.py:

#!/usr/bin/env python
import psyco; psyco.full()
from math import sqrt, ceil
import numpy as np

def rwh_primes(n):
    # /programming/2068372/fastest-way-to-list-all-primes-below-n-in-python/3035188#3035188
    """ Returns  a list of primes < n """
    sieve = [True] * n
    for i in xrange(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 xrange(3,n,2) if sieve[i]]

def rwh_primes1(n):
    # /programming/2068372/fastest-way-to-list-all-primes-below-n-in-python/3035188#3035188
    """ Returns  a list of primes < n """
    sieve = [True] * (n/2)
    for i in xrange(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 xrange(1,n/2) if sieve[i]]

def rwh_primes2(n):
    # /programming/2068372/fastest-way-to-list-all-primes-below-n-in-python/3035188#3035188
    """ Input n>=6, Returns a list of primes, 2 <= p < 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 xrange(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 xrange(1,n/3-correction) if sieve[i]]

def sieve_wheel_30(N):
    # http://zerovolt.com/?p=88
    ''' Returns a list of primes <= N using wheel criterion 2*3*5 = 30

Copyright 2009 by zerovolt.com
This code is free for non-commercial purposes, in which case you can just leave this comment as a credit for my work.
If you need this code for commercial purposes, please contact me by sending an email to: info [at] zerovolt [dot] com.'''
    __smallp = ( 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59,
    61, 67, 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137, 139,
    149, 151, 157, 163, 167, 173, 179, 181, 191, 193, 197, 199, 211, 223, 227,
    229, 233, 239, 241, 251, 257, 263, 269, 271, 277, 281, 283, 293, 307, 311,
    313, 317, 331, 337, 347, 349, 353, 359, 367, 373, 379, 383, 389, 397, 401,
    409, 419, 421, 431, 433, 439, 443, 449, 457, 461, 463, 467, 479, 487, 491,
    499, 503, 509, 521, 523, 541, 547, 557, 563, 569, 571, 577, 587, 593, 599,
    601, 607, 613, 617, 619, 631, 641, 643, 647, 653, 659, 661, 673, 677, 683,
    691, 701, 709, 719, 727, 733, 739, 743, 751, 757, 761, 769, 773, 787, 797,
    809, 811, 821, 823, 827, 829, 839, 853, 857, 859, 863, 877, 881, 883, 887,
    907, 911, 919, 929, 937, 941, 947, 953, 967, 971, 977, 983, 991, 997)

    wheel = (2, 3, 5)
    const = 30
    if N < 2:
        return []
    if N <= const:
        pos = 0
        while __smallp[pos] <= N:
            pos += 1
        return list(__smallp[:pos])
    # make the offsets list
    offsets = (7, 11, 13, 17, 19, 23, 29, 1)
    # prepare the list
    p = [2, 3, 5]
    dim = 2 + N // const
    tk1  = [True] * dim
    tk7  = [True] * dim
    tk11 = [True] * dim
    tk13 = [True] * dim
    tk17 = [True] * dim
    tk19 = [True] * dim
    tk23 = [True] * dim
    tk29 = [True] * dim
    tk1[0] = False
    # help dictionary d
    # d[a , b] = c  ==> if I want to find the smallest useful multiple of (30*pos)+a
    # on tkc, then I need the index given by the product of [(30*pos)+a][(30*pos)+b]
    # in general. If b < a, I need [(30*pos)+a][(30*(pos+1))+b]
    d = {}
    for x in offsets:
        for y in offsets:
            res = (x*y) % const
            if res in offsets:
                d[(x, res)] = y
    # another help dictionary: gives tkx calling tmptk[x]
    tmptk = {1:tk1, 7:tk7, 11:tk11, 13:tk13, 17:tk17, 19:tk19, 23:tk23, 29:tk29}
    pos, prime, lastadded, stop = 0, 0, 0, int(ceil(sqrt(N)))
    # inner functions definition
    def del_mult(tk, start, step):
        for k in xrange(start, len(tk), step):
            tk[k] = False
    # end of inner functions definition
    cpos = const * pos
    while prime < stop:
        # 30k + 7
        if tk7[pos]:
            prime = cpos + 7
            p.append(prime)
            lastadded = 7
            for off in offsets:
                tmp = d[(7, off)]
                start = (pos + prime) if off == 7 else (prime * (const * (pos + 1 if tmp < 7 else 0) + tmp) )//const
                del_mult(tmptk[off], start, prime)
        # 30k + 11
        if tk11[pos]:
            prime = cpos + 11
            p.append(prime)
            lastadded = 11
            for off in offsets:
                tmp = d[(11, off)]
                start = (pos + prime) if off == 11 else (prime * (const * (pos + 1 if tmp < 11 else 0) + tmp) )//const
                del_mult(tmptk[off], start, prime)
        # 30k + 13
        if tk13[pos]:
            prime = cpos + 13
            p.append(prime)
            lastadded = 13
            for off in offsets:
                tmp = d[(13, off)]
                start = (pos + prime) if off == 13 else (prime * (const * (pos + 1 if tmp < 13 else 0) + tmp) )//const
                del_mult(tmptk[off], start, prime)
        # 30k + 17
        if tk17[pos]:
            prime = cpos + 17
            p.append(prime)
            lastadded = 17
            for off in offsets:
                tmp = d[(17, off)]
                start = (pos + prime) if off == 17 else (prime * (const * (pos + 1 if tmp < 17 else 0) + tmp) )//const
                del_mult(tmptk[off], start, prime)
        # 30k + 19
        if tk19[pos]:
            prime = cpos + 19
            p.append(prime)
            lastadded = 19
            for off in offsets:
                tmp = d[(19, off)]
                start = (pos + prime) if off == 19 else (prime * (const * (pos + 1 if tmp < 19 else 0) + tmp) )//const
                del_mult(tmptk[off], start, prime)
        # 30k + 23
        if tk23[pos]:
            prime = cpos + 23
            p.append(prime)
            lastadded = 23
            for off in offsets:
                tmp = d[(23, off)]
                start = (pos + prime) if off == 23 else (prime * (const * (pos + 1 if tmp < 23 else 0) + tmp) )//const
                del_mult(tmptk[off], start, prime)
        # 30k + 29
        if tk29[pos]:
            prime = cpos + 29
            p.append(prime)
            lastadded = 29
            for off in offsets:
                tmp = d[(29, off)]
                start = (pos + prime) if off == 29 else (prime * (const * (pos + 1 if tmp < 29 else 0) + tmp) )//const
                del_mult(tmptk[off], start, prime)
        # now we go back to top tk1, so we need to increase pos by 1
        pos += 1
        cpos = const * pos
        # 30k + 1
        if tk1[pos]:
            prime = cpos + 1
            p.append(prime)
            lastadded = 1
            for off in offsets:
                tmp = d[(1, off)]
                start = (pos + prime) if off == 1 else (prime * (const * pos + tmp) )//const
                del_mult(tmptk[off], start, prime)
    # time to add remaining primes
    # if lastadded == 1, remove last element and start adding them from tk1
    # this way we don't need an "if" within the last while
    if lastadded == 1:
        p.pop()
    # now complete for every other possible prime
    while pos < len(tk1):
        cpos = const * pos
        if tk1[pos]: p.append(cpos + 1)
        if tk7[pos]: p.append(cpos + 7)
        if tk11[pos]: p.append(cpos + 11)
        if tk13[pos]: p.append(cpos + 13)
        if tk17[pos]: p.append(cpos + 17)
        if tk19[pos]: p.append(cpos + 19)
        if tk23[pos]: p.append(cpos + 23)
        if tk29[pos]: p.append(cpos + 29)
        pos += 1
    # remove exceeding if present
    pos = len(p) - 1
    while p[pos] > N:
        pos -= 1
    if pos < len(p) - 1:
        del p[pos+1:]
    # return p list
    return p

def sieveOfEratosthenes(n):
    """sieveOfEratosthenes(n): return the list of the primes < n."""
    # Code from: <dickinsm@gmail.com>, Nov 30 2006
    # http://groups.google.com/group/comp.lang.python/msg/f1f10ced88c68c2d
    if n <= 2:
        return []
    sieve = 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 <end
    using the Sieve of Atkin."""
    # Code by Steve Krenzel, <Sgk284@gmail.com>, 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 xrange(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 xrange((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 xrange(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 xrange((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 xrange(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 xrange(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 xrange(5 >> 1, (int(sqrt(end))+1) >> 1):
        if sieve[n]:
            primes.append((n << 1) + 1)
            aux = (n << 1) + 1
            aux *= aux
            for k in xrange(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 xrange(s, end, 2) if sieve[i >> 1]])

    return primes

def ambi_sieve_plain(n):
    s = range(3, n, 2)
    for m in xrange(3, int(n**0.5)+1, 2): 
        if s[(m-3)/2]: 
            for t in xrange((m*m-3)/2,(n>>1)-1,m):
                s[t]=0
    return [2]+[t for t in s if t>0]

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

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

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

################################################################################
# Using Numpy:
def ambi_sieve(n):
    # http://tommih.blogspot.com/2009/04/fast-prime-number-generator.html
    s = np.arange(3, n, 2)
    for m in xrange(3, int(n ** 0.5)+1, 2): 
        if s[(m-3)/2]: 
            s[(m*m-3)/2::m]=0
    return np.r_[2, s[s>0]]

def primesfrom3to(n):
    # /programming/2068372/fastest-way-to-list-all-primes-below-n-in-python/3035188#3035188
    """ Returns a array of primes, p < n """
    assert n>=2
    sieve = np.ones(n/2, dtype=np.bool)
    for i in xrange(3,int(n**0.5)+1,2):
        if sieve[i/2]:
            sieve[i*i/2::i] = False
    return np.r_[2, 2*np.nonzero(sieve)[0][1::]+1]    

def primesfrom2to(n):
    # /programming/2068372/fastest-way-to-list-all-primes-below-n-in-python/3035188#3035188
    """ Input n>=6, Returns a array of primes, 2 <= p < n """
    sieve = np.ones(n/3 + (n%6==2), dtype=np.bool)
    sieve[0] = False
    for i in xrange(int(n**0.5)/3+1):
        if sieve[i]:
            k=3*i+1|1
            sieve[      ((k*k)/3)      ::2*k] = False
            sieve[(k*k+4*k-2*k*(i&1))/3::2*k] = False
    return np.r_[2,3,((3*np.nonzero(sieve)[0]+1)|1)]

if __name__=='__main__':
    import itertools
    import sys

    def test(f1,f2,num):
        print('Testing {f1} and {f2} return same results'.format(
            f1=f1.func_name,
            f2=f2.func_name))
        if not all([a==b for a,b in itertools.izip_longest(f1(num),f2(num))]):
            sys.exit("Error: %s(%s) != %s(%s)"%(f1.func_name,num,f2.func_name,num))

    n=1000000
    test(sieveOfAtkin,sieveOfEratosthenes,n)
    test(sieveOfAtkin,ambi_sieve,n)
    test(sieveOfAtkin,ambi_sieve_plain,n) 
    test(sieveOfAtkin,sundaram3,n)
    test(sieveOfAtkin,sieve_wheel_30,n)
    test(sieveOfAtkin,primesfrom3to,n)
    test(sieveOfAtkin,primesfrom2to,n)
    test(sieveOfAtkin,rwh_primes,n)
    test(sieveOfAtkin,rwh_primes1,n)         
    test(sieveOfAtkin,rwh_primes2,n)

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

unutbu
источник
4
Если вас интересует не чистый Python-код, вам следует проверить gmpy- у него довольно хорошая поддержка простых чисел с помощью next_primeметода его mpzтипа.
Алекс Мартелли
1
Если вы используете pypy, эти тесты (психологические) кажутся довольно неправильными. Удивительно, но я нашел sieveOfEratosthenes и ambi_sieve_plain самыми быстрыми с pypy. Это то, что я нашел для не пестрых gist.github.com/5bf466bb1ee9e5726a52
Ehsan Kia
1
Если кому-то интересно, как здесь работают функции против PG7.8 в Wikibooks для чистого python без psyco или pypy: для n = 1000000: PG7.8: 4.93 с на цикл; rwh_primes1: 69 мс на цикл; rwh_primes2: 57,1 мс за цикл
потрясающе
8
Можете ли вы обновить это с помощью PyPy, теперь, когда psyco мертв и PyPy заменил его?
noɥʇʎԀʎzɐɹƆ
3
Было бы здорово, если бы эти функции и тайминги могли быть обновлены для python3.
cs95
135

Быстрее и больше памяти чистый код Python:

def primes(n):
    """ Returns  a list of primes < n """
    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 primes1(n):
    """ Returns  a list of primes < n """
    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]]

Быстрее и больше памяти кода:

import numpy
def primesfrom3to(n):
    """ Returns a array of primes, 3 <= p < n """
    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

более быстрое изменение, начинающееся с трети сита:

import numpy
def primesfrom2to(n):
    """ Input n>=6, Returns a array of primes, 2 <= p < n """
    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)]

(Трудно кодируемая) версия чистого кода выше для Python будет выглядеть так:

def primes2(n):
    """ Input n>=6, Returns a list of primes, 2 <= p < n """
    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]]

К сожалению, pure-python не использует более простой и быстрый способ выполнения присваивания, и вызов len()внутри цикла, например, [False]*len(sieve[((k*k)//3)::2*k])слишком медленный. Поэтому мне пришлось импровизировать, чтобы исправить ввод (и избегать больше математики) и сделать некоторую экстремальную (и болезненную) математику.

Лично я думаю, что это позор, что numpy (который так широко используется) не является частью стандартной библиотеки Python, и что улучшения в синтаксисе и скорости, кажется, полностью игнорируются разработчиками Python.

Роберт Уильям Хэнкс
источник
2
Numpy теперь совместим с Python 3. Тот факт, что его нет в стандартной библиотеке, хорош тем, что у них может быть свой цикл выпуска.
Адам
чтобы просто хранить двоичные значения в массиве, я предлагаю bitarray- как здесь используется (для простейшего простого решета; здесь нет соперника в гонке!) stackoverflow.com/questions/31120986/…
hiro главный герой
При приведении в primesfrom2to()методе, должно ли деление быть внутри скобок?
355дней113
3
Для полной версии Python, совместимой с Python 3, перейдите по этой ссылке: stackoverflow.com/a/33356284/2482582
Мобиус
1
Священные окурки это присоска быстро.
Скотт
42

Там довольно аккуратный пример из Python Cookbook здесь - самый быстрый вариант , предложенный по этому URL является:

import itertools
def erat2( ):
    D = {  }
    yield 2
    for q in itertools.islice(itertools.count(3), 0, None, 2):
        p = D.pop(q, None)
        if p is None:
            D[q*q] = q
            yield q
        else:
            x = p + q
            while x in D or not (x&1):
                x += p
            D[x] = p

так что дало бы

def get_primes_erat(n):
  return list(itertools.takewhile(lambda p: p<n, erat2()))

Измеряя в приглашении оболочки (как я предпочитаю делать) с этим кодом в pri.py, я наблюдаю:

$ python2.5 -mtimeit -s'import pri' 'pri.get_primes(1000000)'
10 loops, best of 3: 1.69 sec per loop
$ python2.5 -mtimeit -s'import pri' 'pri.get_primes_erat(1000000)'
10 loops, best of 3: 673 msec per loop

Похоже, что решение Cookbook более чем в два раза быстрее.

Алекс Мартелли
источник
1
@jbochi, пожалуйста - но посмотрите на этот URL, включая кредиты: нам потребовалось десять человек , чтобы все вместе усовершенствовать код до этого момента, включая таких знаменитостей Python, как Tim Peters и Raymond Hettinger (я написал окончательный текст рецепта с тех пор, как я отредактировал печатную поваренную книгу, но с точки зрения кодирования мой вклад был на уровне других) - в конце концов, это действительно тонкий и точно настроенный код, и это не удивительно! -)
Алекс Мартелли
@ Алекс: Зная, что твой код "всего лишь" в два раза быстрее моего, я очень горжусь этим. :) URL тоже было очень интересно читать. Еще раз спасибо.
Jbochi
И это можно сделать еще быстрее с изменением незначительной: см stackoverflow.com/questions/2211990/...
tzot
1
... И это можно сделать еще быстрее с дополнительным ускорением ~ 1,2x-1,3x, резким сокращением объема памяти с O (n) до O (sqrt (n)) и улучшением эмпирической временной сложности, отложив добавление простые числа до тех пор, пока их квадрат не виден на входе. Проверьте это здесь .
Уилл Несс
28

Используя Сито Сундарама , я думаю, что побил рекорд чистого Python:

def sundaram3(max_n):
    numbers = range(3, max_n+1, 2)
    half = (max_n)//2
    initial = 4

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

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

Comparasion:

C:\USERS>python -m timeit -n10 -s "import get_primes" "get_primes.get_primes_erat(1000000)"
10 loops, best of 3: 710 msec per loop

C:\USERS>python -m timeit -n10 -s "import get_primes" "get_primes.daniel_sieve_2(1000000)"
10 loops, best of 3: 435 msec per loop

C:\USERS>python -m timeit -n10 -s "import get_primes" "get_primes.sundaram3(1000000)"
10 loops, best of 3: 327 msec per loop
jbochi
источник
1
Мне удалось ускорить вашу функцию примерно на 20%, добавив «zero = 0» вверху функции, а затем заменив лямбду в вашем фильтре на «zero .__ sub__». Не самый красивый код в мире, но немного быстрее :)
Truppo
1
@ truppo: Спасибо за ваш комментарий! Я только что понял, что передача Noneвместо оригинальной функции работает, и это даже быстрее, чемzero.__sub__
jbochi
7
Знаете ли вы, что если вы сдадите, sundaram3(9)он вернется [2, 3, 5, 7, 9]? Это , кажется, делает это с многочисленным - возможно , все - нечетными числами (даже если они не являются простыми)
wrhall
1
у него есть проблема: sundaram3 (7071) включает 7071, хотя он не прост
bigOther
18

Алгоритм быстрый, но у него есть серьезный недостаток:

>>> sorted(get_primes(530))
[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73,
79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137, 139, 149, 151, 157, 163,
167, 173, 179, 181, 191, 193, 197, 199, 211, 223, 227, 229, 233, 239, 241, 251,
257, 263, 269, 271, 277, 281, 283, 293, 307, 311, 313, 317, 331, 337, 347, 349,
353, 359, 367, 373, 379, 383, 389, 397, 401, 409, 419, 421, 431, 433, 439, 443,
449, 457, 461, 463, 467, 479, 487, 491, 499, 503, 509, 521, 523, 527, 529]
>>> 17*31
527
>>> 23*23
529

Вы предполагаете, numbers.pop()что вернет наименьшее число в наборе, но это совсем не гарантировано. Наборы неупорядочены, pop()удаляют и возвращают произвольный элемент, поэтому его нельзя использовать для выбора следующего простого числа из оставшихся чисел.

STH
источник
17

Для действительно быстрого решения с достаточно большим N было бы загрузить предварительно вычисленный список простых чисел , сохранить его как кортеж и сделать что-то вроде:

for pos,i in enumerate(primes):
    if i > N:
        print primes[:pos]

Если бы N > primes[-1] только тогда вычислите больше простых чисел и сохраните новый список в своем коде, так что в следующий раз он будет таким же быстрым.

Всегда мыслить нестандартно.

Kimvais
источник
9
Чтобы быть справедливым, однако, вы должны были бы посчитать время загрузки, распаковки и форматирования простых чисел и сравнить их со временем для генерации простых чисел с использованием алгоритма - любой из этих алгоритмов может легко записать результаты в файл для последующего использования. использовать. Я думаю, что в этом случае, при наличии достаточного количества памяти для вычисления всех простых чисел, меньших 982 451 653, решение по крошке все равно будет быстрее.
Даниэль Г
3
@ Даниэль правильно. Однако магазин того, что у вас есть, и продолжайте, когда это необходимо, все еще стоит ...
Кимвайс
@ Даниель Г. Я думаю, время загрузки не имеет значения. Разве это не правда о генерации чисел, поэтому вам следует принять во внимание алгоритм, используемый для создания загружаемого списка. И любая временная сложность будет игнорировать однократную передачу файла, учитывая его O (n).
Росс
В разделе часто задаваемых вопросов для главной страницы UTM предлагается вычислять небольшие простые числа быстрее, чем считывать их с диска (вопрос в том, что означает небольшое число).
Бэтмен
12

Если вы не хотите изобретать велосипед, вы можете установить библиотеку символических математических файлов sympy (да, она совместима с Python 3)

pip install sympy

И использовать функцию простого ряда

from sympy import sieve
primes = list(sieve.primerange(1, 10**6))
Полковник паника
источник
9

Если вы принимаете itertools, но не NumPy, вот адаптация rwh_primes2 для Python 3, которая работает примерно в два раза быстрее на моей машине. Единственным существенным изменением является использование байтового массива вместо списка для логического значения и использование сжатия вместо понимания списка для построения окончательного списка. (Я бы добавил это как комментарий, как moarningsun, если бы смог.)

import itertools
izip = itertools.zip_longest
chain = itertools.chain.from_iterable
compress = itertools.compress
def rwh_primes2_python3(n):
    """ Input n>=6, Returns a list of primes, 2 <= p < n """
    zero = bytearray([False])
    size = n//3 + (n % 6 == 2)
    sieve = bytearray([True]) * size
    sieve[0] = False
    for i in range(int(n**0.5)//3+1):
      if sieve[i]:
        k=3*i+1|1
        start = (k*k+4*k-2*k*(i&1))//3
        sieve[(k*k)//3::2*k]=zero*((size - (k*k)//3 - 1) // (2 * k) + 1)
        sieve[  start ::2*k]=zero*((size -   start  - 1) // (2 * k) + 1)
    ans = [2,3]
    poss = chain(izip(*[range(i, n, 6) for i in (1,5)]))
    ans.extend(compress(poss, sieve))
    return ans

Сравнения:

>>> timeit.timeit('primes.rwh_primes2(10**6)', setup='import primes', number=1)
0.0652179726976101
>>> timeit.timeit('primes.rwh_primes2_python3(10**6)', setup='import primes', number=1)
0.03267321276325674

а также

>>> timeit.timeit('primes.rwh_primes2(10**8)', setup='import primes', number=1)
6.394284538007014
>>> timeit.timeit('primes.rwh_primes2_python3(10**8)', setup='import primes', number=1)
3.833829450302801
Джейсон
источник
Очень крутая реализация. :)
Криш
7

Поучительно написать свой собственный основной код поиска, но также полезно иметь быструю надежную библиотеку под рукой. Я написал обертку вокруг primesieve библиотеки C ++ , назвал ее primesieve-python

Попробуй это pip install primesieve

import primesieve
primes = primesieve.generate_primes(10**8)

Мне было бы любопытно увидеть скорость по сравнению.

Полковник паника
источник
Это не совсем то, что заказал OP, но я не понимаю, почему понизили. Это решение за 2,8 секунды в отличие от других внешних модулей. Я заметил в источнике, что он резьбовой, есть тесты на то, насколько хорошо он масштабируется?
ljetibo
@ljetibo ура. Узкое место, кажется, копирует вектор C ++ в список Python, поэтому count_primesфункция намного быстрее, чемgenerate_primes
полковник Паник
На моем компьютере он может сгенерировать простые числа до 1e8 (это дает MemoryError для 1e9) и считать простые числа до 1e10. @HappyLeapSecond выше сравнивает алгоритмы для 1e6
полковник Panic
7

Вот две обновленные (чистый Python 3.6) версии одной из самых быстрых функций,

from itertools import compress

def rwh_primes1v1(n):
    """ Returns  a list of primes < n for n > 2 """
    sieve = bytearray([True]) * (n//2)
    for i in range(3,int(n**0.5)+1,2):
        if sieve[i//2]:
            sieve[i*i//2::i] = bytearray((n-i*i-1)//(2*i)+1)
    return [2,*compress(range(3,n,2), sieve[1:])]

def rwh_primes1v2(n):
    """ Returns a list of primes < n for n > 2 """
    sieve = bytearray([True]) * (n//2+1)
    for i in range(1,int(n**0.5)//2+1):
        if sieve[i]:
            sieve[2*i*(i+1)::2*i+1] = bytearray((n//2-2*i*(i+1))//(2*i+1)+1)
    return [2,*compress(range(3,n,2), sieve[1:])]
Бруно Астролино
источник
1
В Python 3 я использовал эту функцию stackoverflow.com/a/3035188/7799269, но заменил / на // и xrange с дальностью, и они казались намного быстрее, чем эта.
Samerivertwice
4

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

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 <9,080,191 для a = 2,3,37, и 73 достаточно, чтобы решить, является ли N составным или нет.

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

Руджеро Спирман
источник
1
Спасибо за тест простоты Миллера-Рабина, но этот код на самом деле медленнее и не дает правильных результатов. 37 простое и не проходит тест.
Jbochi
Я думаю, что 37 является одним из особых случаев, мой плохой. Я надеялся на детерминированную версию, хотя :)
Ruggiero Spearman
Для Рабина Миллера нет особого случая.
Бессмысленно
2
Вы неправильно прочитали статью. Это 31, а не 37. Вот почему ваша реализация терпит неудачу.
Логан
4

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

Дейв В. Смит
источник
3
Или загрузите их отсюда primes.utm.edu/lists/small/millions , но идея состоит в том, чтобы протестировать ограничение python и посмотреть, появляется ли красивый код в результате оптимизации.
Jbochi
4

Вот код, который я обычно использую для генерации простых чисел в 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)

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

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

MAK
источник
3

Для самого быстрого кода, решение NumPy является лучшим. Однако по чисто академическим причинам я публикую свою версию на чистом 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

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

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

импорт математики
импорт NumPy
def prime6 (до):
    простые числа = numpy.arange (3, ДО + 1,2)
    IsPrime = numpy.ones ((до-1) / 2, DTYPE = BOOL)
    для фактора в простых числах [: int (math.sqrt (upto))]:
        if isprime [(factor-2) / 2]: isprime [(factor * 3-2) / 2: (upto-1) / 2: factor] = 0
    вернуть numpy.insert (простые [isprime], 0,2)

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

nolfonzo
источник
upto=10**6: primesfrom2to()- 7 мс; prime6()- 12 мс ideone.com/oDg2Y
JFS
3

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

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 мс для лучшего решения. Сравните это !!!

lifolofi
источник
3

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

Победителями являются функции, которые используют библиотеку NumPy,

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

Результат расчета времени

Образец кода

Полный код в моем репозитории 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):
    # /programming/579310/formatting-long-numbers-as-strings-in-python?answertab=active#tab-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()
Бруно Адели
источник
2
для сравнения алгоритмических характеристик лучше построить график в масштабе log-log .
Уилл Несс
3

Для 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]]
SmartManoj
источник
3

Самое быстрое простое сито в Pure Python :

from itertools import compress

def half_sieve(n):
    """
    Returns a list of prime numbers less than `n`.
    """
    if n <= 2:
        return []
    sieve = bytearray([True]) * (n // 2)
    for i in range(3, int(n ** 0.5) + 1, 2):
        if sieve[i // 2]:
            sieve[i * i // 2::i] = bytearray((n - i * i - 1) // (2 * i) + 1)
    primes = list(compress(range(1, n, 2), sieve))
    primes[0] = 2
    return primes

Я оптимизировал Sieve of Eratosthenes для скорости и памяти.

эталонный тест

from time import clock
import platform

def benchmark(iterations, limit):
    start = clock()
    for x in range(iterations):
        half_sieve(limit)
    end = clock() - start
    print(f'{end/iterations:.4f} seconds for primes < {limit}')

if __name__ == '__main__':
    print(platform.python_version())
    print(platform.platform())
    print(platform.processor())
    it = 10
    for pw in range(4, 9):
        benchmark(it, 10**pw)

Вывод

>>> 3.6.7
>>> Windows-10-10.0.17763-SP0
>>> Intel64 Family 6 Model 78 Stepping 3, GenuineIntel
>>> 0.0003 seconds for primes < 10000
>>> 0.0021 seconds for primes < 100000
>>> 0.0204 seconds for primes < 1000000
>>> 0.2389 seconds for primes < 10000000
>>> 2.6702 seconds for primes < 100000000
MrSeeker
источник
2

Впервые использую 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 секунд!

smac89
источник
2

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

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

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_pri (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_pri (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_pri (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
АБРИ
источник
чтобы визуализировать результаты своего теста, нанесите их на график в логарифмическом масштабе, чтобы увидеть и сравнить эмпирические порядки роста .
Уилл Несс
@ Будет спасибо за вклад, я буду иметь это в виду в следующий раз, когда мне понадобится такое сравнение
ABri
1

Вот небольшая версия Sieve of Eratosthenes, имеющая как хорошую сложность (ниже, чем сортировка массива длины n), так и векторизацию. По сравнению с @unutbu раз это так же быстро, как пакеты с 46 микросекундами, чтобы найти все простые числа меньше миллиона.

import numpy as np 
def generate_primes(n):
    is_prime = np.ones(n+1,dtype=bool)
    is_prime[0:2] = False
    for i in range(int(n**0.5)+1):
        if is_prime[i]:
            is_prime[i**2::i]=False
    return np.where(is_prime)[0]

Тайминги:

import time    
for i in range(2,10):
    timer =time.time()
    generate_primes(10**i)
    print('n = 10^',i,' time =', round(time.time()-timer,6))

>> n = 10^ 2  time = 5.6e-05
>> n = 10^ 3  time = 6.4e-05
>> n = 10^ 4  time = 0.000114
>> n = 10^ 5  time = 0.000593
>> n = 10^ 6  time = 0.00467
>> n = 10^ 7  time = 0.177758
>> n = 10^ 8  time = 1.701312
>> n = 10^ 9  time = 19.322478
Питер Молгаард Паллесен
источник
1

Я обновил большую часть кода для Python 3 и бросил его на perfplot ( мой проект), чтобы увидеть, какой из них действительно самый быстрый. Оказывается, что, по большому счету n, primesfrom{2,3}toвозьмите торт:

введите описание изображения здесь


Код для воспроизведения сюжета:

import perfplot
from math import sqrt, ceil
import numpy as np
import sympy


def rwh_primes(n):
    # /programming/2068372/fastest-way-to-list-all-primes-below-n-in-python/3035188#3035188
    """ Returns  a list of primes < n """
    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):
    # /programming/2068372/fastest-way-to-list-all-primes-below-n-in-python/3035188#3035188
    """ Returns  a list of primes < n """
    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):
    # /programming/2068372/fastest-way-to-list-all-primes-below-n-in-python/3035188#3035188
    """Input n>=6, Returns a list of primes, 2 <= p < n"""
    assert n >= 6
    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]]


def sieve_wheel_30(N):
    # http://zerovolt.com/?p=88
    """ Returns a list of primes <= N using wheel criterion 2*3*5 = 30

Copyright 2009 by zerovolt.com
This code is free for non-commercial purposes, in which case you can just leave this comment as a credit for my work.
If you need this code for commercial purposes, please contact me by sending an email to: info [at] zerovolt [dot] com."""
    __smallp = (
        2,
        3,
        5,
        7,
        11,
        13,
        17,
        19,
        23,
        29,
        31,
        37,
        41,
        43,
        47,
        53,
        59,
        61,
        67,
        71,
        73,
        79,
        83,
        89,
        97,
        101,
        103,
        107,
        109,
        113,
        127,
        131,
        137,
        139,
        149,
        151,
        157,
        163,
        167,
        173,
        179,
        181,
        191,
        193,
        197,
        199,
        211,
        223,
        227,
        229,
        233,
        239,
        241,
        251,
        257,
        263,
        269,
        271,
        277,
        281,
        283,
        293,
        307,
        311,
        313,
        317,
        331,
        337,
        347,
        349,
        353,
        359,
        367,
        373,
        379,
        383,
        389,
        397,
        401,
        409,
        419,
        421,
        431,
        433,
        439,
        443,
        449,
        457,
        461,
        463,
        467,
        479,
        487,
        491,
        499,
        503,
        509,
        521,
        523,
        541,
        547,
        557,
        563,
        569,
        571,
        577,
        587,
        593,
        599,
        601,
        607,
        613,
        617,
        619,
        631,
        641,
        643,
        647,
        653,
        659,
        661,
        673,
        677,
        683,
        691,
        701,
        709,
        719,
        727,
        733,
        739,
        743,
        751,
        757,
        761,
        769,
        773,
        787,
        797,
        809,
        811,
        821,
        823,
        827,
        829,
        839,
        853,
        857,
        859,
        863,
        877,
        881,
        883,
        887,
        907,
        911,
        919,
        929,
        937,
        941,
        947,
        953,
        967,
        971,
        977,
        983,
        991,
        997,
    )
    # wheel = (2, 3, 5)
    const = 30
    if N < 2:
        return []
    if N <= const:
        pos = 0
        while __smallp[pos] <= N:
            pos += 1
        return list(__smallp[:pos])
    # make the offsets list
    offsets = (7, 11, 13, 17, 19, 23, 29, 1)
    # prepare the list
    p = [2, 3, 5]
    dim = 2 + N // const
    tk1 = [True] * dim
    tk7 = [True] * dim
    tk11 = [True] * dim
    tk13 = [True] * dim
    tk17 = [True] * dim
    tk19 = [True] * dim
    tk23 = [True] * dim
    tk29 = [True] * dim
    tk1[0] = False
    # help dictionary d
    # d[a , b] = c  ==> if I want to find the smallest useful multiple of (30*pos)+a
    # on tkc, then I need the index given by the product of [(30*pos)+a][(30*pos)+b]
    # in general. If b < a, I need [(30*pos)+a][(30*(pos+1))+b]
    d = {}
    for x in offsets:
        for y in offsets:
            res = (x * y) % const
            if res in offsets:
                d[(x, res)] = y
    # another help dictionary: gives tkx calling tmptk[x]
    tmptk = {1: tk1, 7: tk7, 11: tk11, 13: tk13, 17: tk17, 19: tk19, 23: tk23, 29: tk29}
    pos, prime, lastadded, stop = 0, 0, 0, int(ceil(sqrt(N)))

    # inner functions definition
    def del_mult(tk, start, step):
        for k in range(start, len(tk), step):
            tk[k] = False

    # end of inner functions definition
    cpos = const * pos
    while prime < stop:
        # 30k + 7
        if tk7[pos]:
            prime = cpos + 7
            p.append(prime)
            lastadded = 7
            for off in offsets:
                tmp = d[(7, off)]
                start = (
                    (pos + prime)
                    if off == 7
                    else (prime * (const * (pos + 1 if tmp < 7 else 0) + tmp)) // const
                )
                del_mult(tmptk[off], start, prime)
        # 30k + 11
        if tk11[pos]:
            prime = cpos + 11
            p.append(prime)
            lastadded = 11
            for off in offsets:
                tmp = d[(11, off)]
                start = (
                    (pos + prime)
                    if off == 11
                    else (prime * (const * (pos + 1 if tmp < 11 else 0) + tmp)) // const
                )
                del_mult(tmptk[off], start, prime)
        # 30k + 13
        if tk13[pos]:
            prime = cpos + 13
            p.append(prime)
            lastadded = 13
            for off in offsets:
                tmp = d[(13, off)]
                start = (
                    (pos + prime)
                    if off == 13
                    else (prime * (const * (pos + 1 if tmp < 13 else 0) + tmp)) // const
                )
                del_mult(tmptk[off], start, prime)
        # 30k + 17
        if tk17[pos]:
            prime = cpos + 17
            p.append(prime)
            lastadded = 17
            for off in offsets:
                tmp = d[(17, off)]
                start = (
                    (pos + prime)
                    if off == 17
                    else (prime * (const * (pos + 1 if tmp < 17 else 0) + tmp)) // const
                )
                del_mult(tmptk[off], start, prime)
        # 30k + 19
        if tk19[pos]:
            prime = cpos + 19
            p.append(prime)
            lastadded = 19
            for off in offsets:
                tmp = d[(19, off)]
                start = (
                    (pos + prime)
                    if off == 19
                    else (prime * (const * (pos + 1 if tmp < 19 else 0) + tmp)) // const
                )
                del_mult(tmptk[off], start, prime)
        # 30k + 23
        if tk23[pos]:
            prime = cpos + 23
            p.append(prime)
            lastadded = 23
            for off in offsets:
                tmp = d[(23, off)]
                start = (
                    (pos + prime)
                    if off == 23
                    else (prime * (const * (pos + 1 if tmp < 23 else 0) + tmp)) // const
                )
                del_mult(tmptk[off], start, prime)
        # 30k + 29
        if tk29[pos]:
            prime = cpos + 29
            p.append(prime)
            lastadded = 29
            for off in offsets:
                tmp = d[(29, off)]
                start = (
                    (pos + prime)
                    if off == 29
                    else (prime * (const * (pos + 1 if tmp < 29 else 0) + tmp)) // const
                )
                del_mult(tmptk[off], start, prime)
        # now we go back to top tk1, so we need to increase pos by 1
        pos += 1
        cpos = const * pos
        # 30k + 1
        if tk1[pos]:
            prime = cpos + 1
            p.append(prime)
            lastadded = 1
            for off in offsets:
                tmp = d[(1, off)]
                start = (
                    (pos + prime)
                    if off == 1
                    else (prime * (const * pos + tmp)) // const
                )
                del_mult(tmptk[off], start, prime)
    # time to add remaining primes
    # if lastadded == 1, remove last element and start adding them from tk1
    # this way we don't need an "if" within the last while
    if lastadded == 1:
        p.pop()
    # now complete for every other possible prime
    while pos < len(tk1):
        cpos = const * pos
        if tk1[pos]:
            p.append(cpos + 1)
        if tk7[pos]:
            p.append(cpos + 7)
        if tk11[pos]:
            p.append(cpos + 11)
        if tk13[pos]:
            p.append(cpos + 13)
        if tk17[pos]:
            p.append(cpos + 17)
        if tk19[pos]:
            p.append(cpos + 19)
        if tk23[pos]:
            p.append(cpos + 23)
        if tk29[pos]:
            p.append(cpos + 29)
        pos += 1
    # remove exceeding if present
    pos = len(p) - 1
    while p[pos] > N:
        pos -= 1
    if pos < len(p) - 1:
        del p[pos + 1 :]
    # return p list
    return p


def sieve_of_eratosthenes(n):
    """sieveOfEratosthenes(n): return the list of the primes < n."""
    # Code from: <dickinsm@gmail.com>, 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 sieve_of_atkin(end):
    """return a list of all the prime numbers <end using the Sieve of Atkin."""
    # Code by Steve Krenzel, <Sgk284@gmail.com>, 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 ambi_sieve_plain(n):
    s = list(range(3, n, 2))
    for m in range(3, int(n ** 0.5) + 1, 2):
        if s[(m - 3) // 2]:
            for t in range((m * m - 3) // 2, (n >> 1) - 1, m):
                s[t] = 0
    return [2] + [t for t in s if t > 0]


def sundaram3(max_n):
    # /programming/2068372/fastest-way-to-list-all-primes-below-n-in-python/2073279#2073279
    numbers = 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] + filter(None, numbers)


# Using Numpy:
def ambi_sieve(n):
    # http://tommih.blogspot.com/2009/04/fast-prime-number-generator.html
    s = np.arange(3, n, 2)
    for m in range(3, int(n ** 0.5) + 1, 2):
        if s[(m - 3) // 2]:
            s[(m * m - 3) // 2::m] = 0
    return np.r_[2, s[s > 0]]


def primesfrom3to(n):
    # /programming/2068372/fastest-way-to-list-all-primes-below-n-in-python/3035188#3035188
    """ Returns an array of primes, p < n """
    assert n >= 2
    sieve = np.ones(n // 2, dtype=np.bool)
    for i in range(3, int(n ** 0.5) + 1, 2):
        if sieve[i // 2]:
            sieve[i * i // 2::i] = False
    return np.r_[2, 2 * np.nonzero(sieve)[0][1::] + 1]


def primesfrom2to(n):
    # /programming/2068372/fastest-way-to-list-all-primes-below-n-in-python/3035188#3035188
    """ Input n>=6, Returns an array of primes, 2 <= p < n """
    assert n >= 6
    sieve = np.ones(n // 3 + (n % 6 == 2), dtype=np.bool)
    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
            sieve[(k * k + 4 * k - 2 * k * (i & 1)) // 3::2 * k] = False
    return np.r_[2, 3, ((3 * np.nonzero(sieve)[0] + 1) | 1)]


def sympy_sieve(n):
    return list(sympy.sieve.primerange(1, n))


perfplot.save(
    "prime.png",
    setup=lambda n: n,
    kernels=[
        rwh_primes,
        rwh_primes1,
        rwh_primes2,
        sieve_wheel_30,
        sieve_of_eratosthenes,
        sieve_of_atkin,
        # ambi_sieve_plain,
        # sundaram3,
        ambi_sieve,
        primesfrom3to,
        primesfrom2to,
        sympy_sieve,
    ],
    n_range=[2 ** k for k in range(3, 25)],
    logx=True,
    logy=True,
    xlabel="n",
)
Нико Шлёмер
источник
0

Я предполагаю, что самый быстрый из всех способов - это жестко закодировать простые числа в вашем коде.

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

Конечно, это работает, только если вы знаете верхнюю границу N во время компиляции, но это относится к (почти) всем проблемам Эйлера проекта.

 

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

akuhn
источник
0

Извините, что беспокою, но erat2 () имеет серьезный недостаток в алгоритме.

При поиске следующей композиции нам нужно проверить только нечетные числа. q, p оба нечетны; тогда q + p четно и не нуждается в проверке, но q + 2 * p всегда нечетно. Это исключает проверку «если даже» в цикле while и экономит около 30% времени выполнения.

Пока мы на этом: вместо элегантного метода «D.pop (q, None)» get и delete используйте «if q в D: p = D [q], del D [q]», что в два раза быстрее ! По крайней мере, на моей машине (P3-1Ghz). Поэтому я предлагаю эту реализацию этого умного алгоритма:

def erat3( ):
    from itertools import islice, count

    # q is the running integer that's checked for primeness.
    # yield 2 and no other even number thereafter
    yield 2
    D = {}
    # no need to mark D[4] as we will test odd numbers only
    for q in islice(count(3),0,None,2):
        if q in D:                  #  is composite
            p = D[q]
            del D[q]
            # q is composite. p=D[q] is the first prime that
            # divides it. Since we've reached q, we no longer
            # need it in the map, but we'll mark the next
            # multiple of its witnesses to prepare for larger
            # numbers.
            x = q + p+p        # next odd(!) multiple
            while x in D:      # skip composites
                x += p+p
            D[x] = p
        else:                  # is prime
            # q is a new prime.
            # Yield it and mark its first multiple that isn't
            # already marked in previous iterations.
            D[q*q] = q
            yield q
user1016274
источник
Относительно отложенного добавления простых чисел в dict (до тех пор, пока во входных данных не будет виден квадрат простого числа), см. stackoverflow.com/a/10733621/849891 .
Будет Несс
0

Самый быстрый метод, который я пробовал до сих пор, основан на функции поваренной книги Pythonerat2 :

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:
            x = q + 2*p
            while x in D:
                x += 2*p
            D[x] = p

Смотрите этот ответ для объяснения ускорения.

tzot
источник
0

Я могу опоздать на вечеринку, но мне придется добавить свой код для этого. Он использует приблизительно n / 2 в пространстве, потому что нам не нужно хранить четные числа, и я также использую модуль python bitarray, дополнительно сокращая потребление памяти и позволяя вычислять все простые числа до 1 000 000 000

from bitarray import bitarray
def primes_to(n):
    size = n//2
    sieve = bitarray(size)
    sieve.setall(1)
    limit = int(n**0.5)
    for i in range(1,limit):
        if sieve[i]:
            val = 2*i+1
            sieve[(i+i*val)::val] = 0
    return [2] + [2*i+1 for i, v in enumerate(sieve) if v and i > 0]

python -m timeit -n10 -s "import euler" "euler.primes_to(1000000000)"
10 loops, best of 3: 46.5 sec per loop

Это было запущено на 64-битной 2.4GHZ MAC OSX 10.8.3

Cobie
источник
1
размещение одного времени для неизвестной машины ничего не говорит. Принятый ответ здесь гласит: «без psyco, для n = 1000000 rwh_primes2 был самым быстрым». Так что, если бы вы предоставили свои сроки для этого кода, а также для своего кода на той же машине, а также в 2, 4, 10 млн., То это было бы гораздо более информативно.
Уилл Несс
-1, этот код зависит от особенностей bitarray, реализованного в C, поэтому код работает быстро, так как большая часть работы выполняется в нативном коде при назначении слайса. BitArray пакет BREAKS стандартного определения для собственных срезов (индексированного по диапазону) для изменяемых последовательностей в том , что она позволяет назначить одну или логическое значение 0/1 True / False для всех элементов среза, в то время как стандартное поведение для чистого Python , кажется, чтобы не допустить этого и разрешить только значение присваивания 0, в этом случае оно рассматривается как деление всех элементов слайса из последовательности / массива.
GordonBGood
продолжение: Если сравнивать вызов нестандартного нативного кода, мы можем также написать пакет генератора последовательностей «fastprimes», основанный на C-коде, например, из простого числа Кима Валиша, и сгенерировать все простые числа из четырех миллиардов плюс 32 -битный диапазон номеров всего за несколько секунд с помощью одного вызова генератора последовательности. При этом также практически не использовалась бы память, поскольку связанный код основан на сегментированном сите Эратосфена и, таким образом, использует только несколько десятков килобайт ОЗУ, и если бы была сгенерирована последовательность, хранилище списков не требовалось.
GordonBGood
0

Со временем я собрал несколько сит с простыми числами. Самый быстрый на моем компьютере это:

from time import time
# 175 ms for all the primes up to the value 10**6
def primes_sieve(limit):
    a = [True] * limit
    a[0] = a[1] = False
    #a[2] = True
    for n in xrange(4, limit, 2):
        a[n] = False
    root_limit = int(limit**.5)+1
    for i in xrange(3,root_limit):
        if a[i]:
            for n in xrange(i*i, limit, 2*i):
                a[n] = False
    return a

LIMIT = 10**6
s=time()
primes = primes_sieve(LIMIT)
print time()-s
Стефан Грюнвальд
источник
0

Я медленно отвечаю на этот вопрос, но это казалось забавным упражнением. Я использую NumPy, который может быть мошенничеством, и я сомневаюсь, что этот метод самый быстрый, но это должно быть ясно. Он фильтрует логический массив, ссылаясь только на его индексы, и извлекает простые числа из индексов всех значений True. Не требуется по модулю.

import numpy as np
def ajs_primes3a(upto):
    mat = np.ones((upto), dtype=bool)
    mat[0] = False
    mat[1] = False
    mat[4::2] = False
    for idx in range(3, int(upto ** 0.5)+1, 2):
        mat[idx*2::idx] = False
    return np.where(mat == True)[0]
Алан Джеймс Салмони
источник
это неверно, например, ajs_primes3a(10)-> array([2, 3, 5, 7, 9]). 9это не простое число
JFS
Вы заметили крайний случай, которого у меня не было - молодец! Проблема была в 'для idx в диапазоне (3, int (до ** 0,5), 2):' который должен быть 'для idx в диапазоне (3, int (до ** 0,5) + 1, 2):'. Спасибо, но теперь это работает.
Алан Джеймс Салмони
Причина была в том, что цикл idx поднялся до «до ** 05», что для случаев до 15. включительно. С 16 и далее он работает нормально. Это был набор крайних случаев, для которых я не проверял. Добавление 1 означает, что оно должно работать для всех чисел.
Алан Джеймс Салмони
Кажется, сейчас работает. Это самое медленное numpyрешение на основе массива, которое возвращает массив. Примечание: истинная реализация Sieve of Eratosthenes не используется по модулю - упоминать об этом не нужно. Вы могли бы использовать mat[idx*idx::idx]вместо mat[idx*2::idx]. И np.nonzero(mat)[0]вместо np.where(mat == True)[0].
JFS
Спасибо JF. Я протестировал против prime6 () и получил результат быстрее до (IIRC) около 250k, когда prime6 () вступил во владение. primesfrom2to () был быстрее. На расстоянии до 20 м ajs_primes3a () занимало 0,034744977951 мс, prime6 () - 0,0222899913788 мс, а primesfrom2to () - 0,0104751586914 мс (та же машина, та же нагрузка, лучшее из 10 таймингов). Честно, лучше, чем я думал.
Алан Джеймс Салмони
0

Вот интересный метод генерации простых чисел (но не самый эффективный) с использованием понимания списка Python:

noprimes = [j for i in range(2, 8) for j in range(i*2, 50, i)]
primes = [x for x in range(2, 50) if x not in noprimes]

Вы можете найти пример и некоторые пояснения прямо здесь

Александр
источник