Najszybszy sposób na wyświetlenie wszystkich liczb pierwszych poniżej N.


357

To najlepszy algorytm, jaki mogłem wymyślić.

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

Czy można to zrobić jeszcze szybciej?

Ten kod ma wadę: ponieważ numbersjest to zestaw nieuporządkowany, nie ma gwarancji, że numbers.pop()usunie z niego najniższą liczbę. Niemniej jednak działa (przynajmniej dla mnie) dla niektórych liczb wejściowych:

>>> 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

Fragment kodu, o którym mowa, jest znacznie szybszy, jeśli liczby zadeklarowane jak liczby = set (zakres (n, 2, -2)). Ale nie można pokonać sundaram3. Dzięki za pytanie.
Shekhar

3
Byłoby miło, gdyby w odpowiedziach mogły znajdować się wersje funkcji w języku Python 3.
Michael Foukarakis,

Na pewno jest do tego biblioteka, więc nie musimy rzucać własnym> xkcd obiecanym Pythonem jest tak prosty jak import antigravity. Czy nie ma czegoś takiego jak require 'prime'; Prime.take(10)(Ruby)?
Pułkownik Panic

2
@ ColonelPanic Tak się składa, że ​​zaktualizowałem github.com/jaredks/pyprimesieve dla Py3 i dodałem do PyPi. Jest to z pewnością szybsze niż te, ale nie rzędy wielkości - więcej niż ~ 5 razy szybsze niż najlepsze wersje numpy.
Jared

3
@ColonelPanic: Wydaje mi się, że edytowanie starych odpowiedzi w celu zauważenia, że ​​się postarzały, jest właściwe, ponieważ czyni je bardziej użytecznym zasobem. Jeśli „zaakceptowana” odpowiedź nie jest już najlepsza, być może edytuj notatkę w pytaniu z aktualizacją z 2015 r., Aby wskazać ludziom obecną najlepszą metodę.
Peter Cordes,

Odpowiedzi:


365

Ostrzeżenie: timeit wyniki mogą się różnić ze względu na różnice w sprzęcie lub wersji Pythona.

Poniżej znajduje się skrypt, który porównuje wiele implementacji:

Wiele dzięki Stephan do wniesienia sieve_wheel_30 moją uwagę. Uznanie dla Roberta Williama Hanksa za primesfrom2to, primesfrom3to, rwh_primes, rwh_primes1 i rwh_primes2.

Spośród testowanych prostych metod Pythona, psyco , dla n = 1000000, rwh_primes1 był najszybciej testowany.

+---------------------+-------+
| 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 |
+---------------------+-------+

Spośród przetestowanych prostych metod Pythona, bez psyco , dla n = 1000000, rwh_primes2 był najszybszy.

+---------------------+-------+
| 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 |
+---------------------+-------+

Ze wszystkich testowanych metod, pozwalających na numpy , dla n = 1000000, primesfrom2to był najszybszym testowanym.

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

Czasy zmierzono za pomocą polecenia:

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

z {method}zastąpionymi przez każdą z nazw metod.

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)

Uruchomienie testów skryptu, aby wszystkie implementacje dały ten sam wynik.


4
Jeśli interesuje Cię nieczytelny kod Pythona, powinieneś sprawdzić gmpy- ma całkiem dobrą obsługę liczb pierwszych, za pomocą next_primemetody tego mpztypu.
Alex Martelli,

1
Jeśli używasz pypy, te testy porównawcze (psyco) wydają się dość odbiegające od normy. Zaskakujące jest to, że sieveOfEratosthenes i ambi_sieve_plain są najszybsze z pypy. Oto, co znalazłem dla tych niepoliczalnych gist.github.com/5bf466bb1ee9e5726a52
Ehsan Kia

1
Jeśli ktoś zastanawia się, jak działają tutaj funkcje w stosunku do PG7.8 Wikibooks dla czystego Pythona bez psyco ani pypy: dla n = 1000000: PG7.8: 4.93 s na pętlę; rwh_primes1: 69 ms na pętlę; rwh_primes2: 57,1 ms na pętlę
okropny

8
Czy możesz to zaktualizować za pomocą PyPy, skoro psyco nie żyje, a PyPy go zastąpił?
noɥʇʎԀʎzɐɹƆ

2
Byłoby wspaniale, gdyby te funkcje i czasy mogły zostać zaktualizowane dla python3.
cs95,

134

Szybszy i bardziej czysty pod względem pamięci czysty kod 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]]

lub zaczynając od półsita

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]]

Szybszy i bardziej logiczny kod numpy:

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

szybsza odmiana zaczynająca się od jednej trzeciej sita:

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)]

Wersja (trudna do kodowania) czysto pythonowego powyższego kodu to:

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]]

Niestety, czysty python nie przyjmuje prostszego i szybszego numerycznego sposobu wykonywania przypisania, a wywoływanie len()wewnątrz pętli, jak w, [False]*len(sieve[((k*k)//3)::2*k])jest zbyt wolne. Musiałem więc improwizować, aby poprawić dane wejściowe (i uniknąć więcej matematyki) i wykonać ekstremalną (i bolesną) matematykę.

Osobiście uważam, że szkoda, że ​​numpy (która jest tak szeroko stosowana) nie jest częścią standardowej biblioteki Pythona, a ulepszenia składni i szybkości wydają się być całkowicie pomijane przez programistów Pythona.


2
Numpy jest teraz kompatybilny z Pythonem 3. Fakt, że nie ma go w standardowej bibliotece, jest dobry, w ten sposób mogą mieć swój własny cykl wydawniczy.
Adam

po prostu przechowuję wartości binarne w tablicy, którą sugeruję bitarray- jak tutaj użyto (dla najprostszego sita głównego; nie ma tu miejsca w wyścigu!) stackoverflow.com/questions/31120986/…
protagonista hiro

Czy podczas odlewania primesfrom2to()metody podział powinien znajdować się w nawiasach?
355durch113

3
Aby uzyskać czystą wersję Pythona zgodną z Pythonem 3, skorzystaj z tego linku: stackoverflow.com/a/33356284/2482582
Moebius

FWIW, wersja pierwszego bloku kodu jest tematem tego pytania . Należy również wspomnieć, że użyłem Python 3 wersje kodu tutaj .
PM 2,

42

Jest to całkiem zgrabny próbka z Python Cookbook tutaj - najszybsza wersja proponowana na tym adresem jest:

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

to by dało

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

Pomiary w wierszu powłoki (jak wolę to zrobić) za pomocą tego kodu w pri.py, obserwuję:

$ 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

więc wygląda na to, że rozwiązanie Cookbook jest ponad dwukrotnie szybsze.


1
@ jbochi, nie ma za co - ale spójrz na ten adres URL, w tym kredyty: dziesięć osób zajęło nam zbiorowe dopracowanie kodu do tego momentu, w tym także wydajnościowe wersje Pythona, takie jak Tim Peters i Raymond Hettinger (napisałem końcowy tekst przepisu, odkąd zredagowałem drukowaną książkę kucharską, ale pod względem kodowania mój wkład był na równi z innymi)) - w końcu jest to naprawdę subtelny i dopracowany kod, i to nie jest zaskakujące! -)
Alex Martelli

@Alex: Wiedząc, że twój kod jest „tylko” dwa razy szybszy niż mój, jestem z tego powodu dumny. :) Adres URL był również bardzo interesujący do przeczytania. Dzięki jeszcze raz.
jbochi

Można to zrobić jeszcze szybciej przy niewielkiej zmianie: patrz stackoverflow.com/questions/2211990/...
tzot

1
... I może być jeszcze szybszy dzięki dodatkowemu przyspieszeniu ~ 1,2x-1,3x, drastycznemu zmniejszeniu powierzchni pamięci z O (n) do O (sqrt (n)) i poprawie złożoności czasu empirycznego poprzez odroczenie dodania wprowadza się do dykt, dopóki ich kwadrat nie będzie widoczny na wejściu. Sprawdź to tutaj .
Czy Ness

28

Używając sita Sundaram , myślę, że pobiłem rekord czysto Pythona:

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)

Porównanie:

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

1
Udało mi się przyspieszyć twoją funkcję o około 20%, dodając „zero = 0” u góry funkcji, a następnie zamieniając lambda w filtrze na „zero .__ sub__”. Nie najładniejszy kod na świecie, ale nieco szybszy :)
truppo

1
@truppo: Dziękujemy za komentarz! Właśnie zdałem sobie sprawę, że przekazywanie Nonezamiast oryginalnej funkcji działa i jest nawet szybsze niżzero.__sub__
jbochi

7
Czy wiesz, że jeśli zdasz sundaram3(9), wróci [2, 3, 5, 7, 9]? Wydaje się, że robi to z licznymi - być może wszystkimi - liczbami nieparzystymi (nawet jeśli nie są liczbami pierwszymi)
wrhall

1
ma problem: sundaram3 (7071) zawiera 7071, ale nie jest liczbą pierwszą
bigOther

18

Algorytm jest szybki, ale ma poważną wadę:

>>> 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

Zakładasz, numbers.pop()że zwróci najmniejszą liczbę w zestawie, ale nie jest to wcale gwarantowane. Zestawy są nieuporządkowane i pop()usuwają i zwracają dowolny element, więc nie można go użyć do wybrania kolejnej liczby pierwszej z pozostałych liczb.


17

Dla naprawdę najszybszego rozwiązania z wystarczająco dużą wartością N byłoby pobranie wstępnie obliczonej listy liczb pierwszych , przechowanie jej jako krotki i zrobienie czegoś takiego:

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

Jeśli N > primes[-1] tylko, to oblicz więcej liczb pierwszych i zapisz nową listę w kodzie, więc następnym razem będzie równie szybka.

Zawsze myśl nieszablonowo.


9
Aby być jednak uczciwym, trzeba liczyć czas pobierania, rozpakowywania i formatowania liczb pierwszych i porównywać to z czasem generowania liczb pierwszych za pomocą algorytmu - każdy z tych algorytmów może łatwo zapisać wyniki do pliku na później posługiwać się. Myślę, że w takim przypadku, biorąc pod uwagę wystarczającą ilość pamięci, aby faktycznie obliczyć wszystkie liczby pierwsze mniejsze niż 982 451 653, rozwiązanie numpy nadal byłoby szybsze.
Daniel G

3
@Daniel poprawnie. Jednak przechowuj to, co masz i kontynuuj, gdy zajdzie taka potrzeba, nadal trwa ...
Kimvais

@Daniel G Myślę, że czas pobierania nie ma znaczenia. Naprawdę nie chodzi o generowanie liczb, więc warto wziąć pod uwagę algorytm użyty do utworzenia listy, którą pobierasz. I za każdym razem złożoność zignorowałaby jeden transfer pliku, biorąc pod uwagę O (n).
Ross

FAQ dla UTM prime stronie sugeruje obliczania małych liczb pierwszych jest szybszy niż czytając je na dysku (kwestia jest to, co małe środki).
Batman

12

Jeśli nie chcesz wymyślać koła od nowa , możesz zainstalować symboliczną sympię biblioteki matematycznej (tak, jest kompatybilna z Python 3)

pip install sympy

I użyj funkcji primerange

from sympy import sieve
primes = list(sieve.primerange(1, 10**6))

8

Jeśli akceptujesz narzędzia itertools, ale nie numpy, oto adaptacja rwh_primes2 dla Python 3, która działa około dwa razy szybciej na moim komputerze. Jedyną istotną zmianą jest użycie bytearray zamiast listy dla wartości logicznej i użycie kompresji zamiast interpretacji listy do zbudowania ostatecznej listy. (Dodałbym to jako komentarz jak moarningsun, gdybym mógł.)

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

Porównania:

>>> 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

i

>>> 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

Bardzo fajne wdrożenie. :)
Krish

7

Pouczające jest pisanie własnego kodu znalezienia głównej, ale przydatne jest również posiadanie szybkiej niezawodnej biblioteki. Napisałem wrapper wokół biblioteki C ++ primesieve , o nazwie primesieve-python

Spróbuj pip install primesieve

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

Byłbym ciekawy, czy prędkość jest porównywana.


Nie jest to dokładnie to, co zamówił OP, ale nie rozumiem, dlaczego głosowanie negatywne. Jest to rozwiązanie 2,8 s, w przeciwieństwie do innych modułów zewnętrznych. Zauważyłem w źródle, że jest wątkowy, mam jakieś testy, jak dobrze się skaluje?
ljetibo

@ljetibo na zdrowie. Wąskim gardłem wydaje się być kopiowanie wektora C ++ na listę Pythona, dlatego count_primesfunkcja jest znacznie szybsza niżgenerate_primes
pułkownik Panic

Na moim komputerze może wygodnie generować liczby pierwsze do 1e8 (daje MemoryError dla 1e9) i liczyć liczby pierwsze do 1e10. @HappyLeapSecond powyżej porównuje algorytmy dla 1e6
pułkownik Panic

7

Oto dwie zaktualizowane (czysto Python 3.6) wersje jednej z najszybszych funkcji,

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
W Pythonie 3 użyłem tej funkcji stackoverflow.com/a/3035188/7799269, ale zastąpiłem / // // i xrange z zakresem i wydawały się znacznie szybsze.
samerivertwice,

4

Deterministyczna implementacja testu pierwszeństwa Millera-Rabina przy założeniu, że 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)

Zgodnie z artykułem na Wikipedii ( http://en.wikipedia.org/wiki/Miller–Rabin_primality_test ) testowanie N <9,080,191 dla a = 2,3,37, a 73 wystarczy, aby zdecydować, czy N jest złożony, czy nie.

I dostosowałem kod źródłowy z probabilistycznej implementacji oryginalnego testu Millera-Rabina znalezionego tutaj: http://en.literateprograms.org/Miller-Rabin_primality_test_(Python)


1
Dziękujemy za test pierwotności Millera-Rabina, ale ten kod jest w rzeczywistości wolniejszy i nie zapewnia poprawnych wyników. 37 jest liczbą pierwszą i nie przechodzi testu.
jbochi

Chyba 37 to jeden ze specjalnych przypadków, mój zły. Miałem jednak nadzieję na wersję deterministyczną :)
Ruggiero Spearman

Nie ma żadnego specjalnego przypadku dla młyna rabinowego.
Misguided

2
Źle przeczytałeś artykuł. To jest 31, a nie 37. Właśnie dlatego twoja implementacja kończy się niepowodzeniem.
Logan,

4

Jeśli masz kontrolę nad N, najszybszym sposobem wyliczenia wszystkich liczb pierwszych jest ich wstępne obliczenie. Poważnie. Wstępne obliczenia to sposób na przeoczenie optymalizacji.


3
Lub pobierz je stąd primes.utm.edu/lists/small/millions , ale pomysł polega na przetestowaniu limitu pythona i sprawdzeniu, czy piękny kod wyłania się z optymalizacji.
jbochi

4

Oto kod, którego zwykle używam do generowania liczb pierwszych w Pythonie:

$ 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)

Nie może konkurować z szybszymi rozwiązaniami zamieszczonymi tutaj, ale przynajmniej jest to czysty python.

Dziękujemy za opublikowanie tego pytania. Naprawdę dużo się dzisiaj nauczyłem.


3

W przypadku najszybszego kodu najlepsze jest rozwiązanie numpy. Jednak ze względów czysto akademickich publikuję moją czystą wersję Pythona, która jest nieco mniej niż 50% szybsza niż wersja książki kucharskiej opublikowana powyżej. Ponieważ tworzę całą listę w pamięci, potrzebujesz wystarczająco dużo miejsca, aby pomieścić wszystko, ale wydaje się, że skaluje się dość dobrze.

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)

A wyniki:

>>>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

Nieco inna implementacja półsita z użyciem Numpy:

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

import matematyki
importuj numpy
def prime6 (upto):
    liczby pierwsze = numpy.arange (3, upto + 1,2)
    isprime = numpy.ones ((upto-1) / 2, dtype = bool)
    dla czynnika w liczbach pierwszych [: int (math.sqrt (upto))]:
        if isprime [(współczynnik-2) / 2]: isprime [(współczynnik * 3-2) / 2: (upto-1) / 2: współczynnik] = 0
    return numpy.insert (primes [isprime], 0,2)

Czy ktoś może to porównać z innymi czasami? Na mojej maszynie wydaje się całkiem porównywalna z innym pół-sitem Numpy.


upto=10**6: primesfrom2to()- 7 ms; prime6()- 12 ms ideone.com/oDg2Y
jfs

3

Wszystko jest napisane i przetestowane. Nie ma więc potrzeby wymyślania nowego koła.

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

daje nam rekord 12,2 milisekund !

10 loops, best of 10: 12.2 msec per loop

Jeśli to nie jest wystarczająco szybkie, możesz wypróbować PyPy:

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

Co skutkuje w:

10 loops, best of 10: 2.03 msec per loop

Odpowiedź z 247 głosami up-up zawiera 15,9 ms dla najlepszego rozwiązania. Porównaj to !!!


3

Przetestowałem niektóre funkcje unutbu , obliczyłem je z liczbą głodnych milionów

Zwycięzcami są funkcje korzystające z biblioteki numpy,

Uwaga : Ciekawe byłoby również wykonanie testu wykorzystania pamięci :)

Wynik czasu obliczeń

Przykładowy kod

Kompletny kod w moim repozytorium 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
aby porównać wydajność algorytmiczną , lepiej wykreślić w skali log-log .
Czy Ness

3

Dla 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]]

3

Najszybsze sito główne w 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

Zoptymalizowałem Sito Eratostenesa pod kątem szybkości i pamięci.

Reper

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)

Wynik

>>> 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

2

Pierwszy raz korzystam z Pythona, więc niektóre metody, które stosuję w tym przypadku, mogą wydawać się nieco kłopotliwe. Właśnie przekonwertowałem mój kod c ++ do Pythona i to właśnie mam (choć odrobinę wolniej w Pythonie)

#!/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

Znaleziono 664579 liczb pierwszych w 12,799119 sekund!

#!/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

Znaleziono 664579 liczb pierwszych w 10,230172 sekund!

#!/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

Znaleziono 664579 liczb pierwszych w 7,113776 sekund!


2

Wiem, że konkurs jest zamknięty na kilka lat. …

Niemniej jednak jest to moja sugestia dotycząca czystego głównego sita pythonowego, opartego na pomijaniu wielokrotności 2, 3 i 5 poprzez zastosowanie odpowiednich kroków podczas przetwarzania sita do przodu. Niemniej jednak jest on wolniejszy dla N <10 ^ 9 niż lepsze rozwiązania Roberta Williama Hanksa rwh_primes2 i rwh_primes1. Dzięki zastosowaniu sita ctypes.c_ushort powyżej 1,5 * 10 ^ 8 w jakiś sposób dostosowuje się do limitów pamięci.

10 ^ 6

$ python -mtimeit -s "import primeSieveSpeedComp" "primeSieveSpeedComp.primeSieveSeq (1000000)" 10 pętli, najlepiej 3: 46,7 ms na pętlę

do porównania: $ python -mtimeit -s "import primeSieveSpeedComp" "primeSieveSpeedComp.rwh_primes1 (1000000)" 10 pętli, najlepiej 3: 43,2 ms na pętlę, aby porównać: $ python -m timeit -s "import primeSieveSpeedComp" "primeSieveSp_primes2" (1000000) „10 pętli, najlepiej 3: 34,5 ms na pętlę

10 ^ 7

$ python -mtimeit -s "import primeSieveSpeedComp" "primeSieveSpeedComp.primeSieveSeq (10000000)" 10 pętli, najlepiej 3: 530 ms na pętlę

do porównania: $ python -mtimeit -s "import primeSieveSpeedComp" "primeSieveSpeedComp.rwh_primes1 (10000000)" 10 pętli, najlepiej 3: 494 ms na pętlę do porównania: $ python -m timeit -s "import primeSieveSpeedComp" "primeSieveSpeedCimesrwh. (10000000) „10 pętli, najlepiej 3: 375 ms na pętlę

10 ^ 8

$ python -mtimeit -s "import primeSieveSpeedComp" "primeSieveSpeedComp.primeSieveSeq (100000000)" 10 pętli, najlepiej 3: 5,55 s na pętlę

do porównania: $ python -mtimeit -s "import primeSieveSpeedComp" "primeSieveSpeedComp.rwh_primes1 (100000000)" 10 pętli, najlepiej 3: 5,33 s na pętlę do porównania: $ python -m timeit -s "import primeSieveSpeedComp" "primeSieveSpeedCimesrwh. (100000000) „10 pętli, najlepiej 3: 3,95 s na pętlę

10 ^ 9

$ python -mtimeit -s "import primeSieveSpeedComp" "primeSieveSpeedComp.primeSieveSeq (1000000000)" 10 pętli, najlepiej 3: 61,2 s na pętlę

do porównania: $ python -mtimeit -n 3 -s "import primeSieveSpeedComp" "primeSieveSpeedComp.rwh_primes1 (1000000000)" 3 pętle, najlepiej 3: 97,8 s na pętlę

do porównania: $ python -m timeit -s "import primeSieveSpeedComp" "primeSieveSpeedComp.rwh_primes2 (1000000000)" 10 pętli, najlepiej 3: 41,9 s na pętlę

Możesz skopiować poniższy kod do ubuntus primeSieveSpeedComp, aby przejrzeć te testy.

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

aby wizualizować wyniki testu, wykreśl je w skali log-log, aby zobaczyć i porównać empiryczne rzędy wzrostu .
Czy Ness

@ Dzięki za wkład, będę o tym pamiętać następnym razem, gdy będę potrzebować takiego porównania
ABri

1

Oto numpy wersja Sita Eratostenesa, mająca zarówno dobrą złożoność (mniejszą niż sortowanie tablicy długości n), jak i wektoryzację. W porównaniu do czasów @unutbu jest to tak szybkie, jak pakiety z 46 mikrosekundami, aby znaleźć wszystkie liczby pierwsze poniżej miliona.

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]

Czasy:

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

Zaktualizowałem wiele kodu Python 3 i rzuciłem go na perfplot (mój projekt), aby zobaczyć, który z nich jest najszybszy. Okazuje się, że za duża n, primesfrom{2,3}towziąć ciasto:

wprowadź opis zdjęcia tutaj


Kod do odtworzenia fabuły:

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

Domyślam się, że najszybszym ze wszystkich sposobów jest kodowanie liczb pierwszych w twoim kodzie.

Dlaczego więc nie napisać powolnego skryptu, który generuje inny plik źródłowy, w którym zapisane są wszystkie liczby, a następnie zaimportować ten plik źródłowy po uruchomieniu rzeczywistego programu.

Oczywiście działa to tylko wtedy, gdy znasz górną granicę N w czasie kompilacji, ale tak jest w przypadku (prawie) wszystkich problemów Eulera projektu.

 

PS: Mogę się mylić, choć iff parsowanie źródła za pomocą stałych liczb stałych jest wolniejsze niż obliczanie ich w pierwszej kolejności, ale o ile wiem, Python działa ze skompilowanych .pycplików, więc czytanie tablicy binarnej ze wszystkimi liczbami pierwszymi do N powinno być krwawe w takim razie szybko.


0

Przepraszam, że przeszkadzam, ale erat2 () ma poważną wadę w algorytmie.

Szukając następnego kompozytu, musimy przetestować tylko liczby nieparzyste. q, p oba są nieparzyste; wtedy q + p jest parzyste i nie musi być testowane, ale q + 2 * p jest zawsze nieparzyste. Eliminuje to test „jeśli nawet” w pętli while i pozwala zaoszczędzić około 30% czasu wykonywania.

Skoro już jesteśmy przy tym: zamiast eleganckiej metody „D.pop (q, None)” pobierz i usuń użyj „jeśli q w D: p = D [q], del D [q]”, która jest dwa razy szybsza ! Przynajmniej na mojej maszynie (P3-1 Ghz). Sugeruję więc implementację tego sprytnego algorytmu:

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

odroczone dodanie liczb pierwszych do dykta (do momentu pojawienia się kwadratu liczby pierwszej na wejściu) patrz stackoverflow.com/a/10733621/849891 .
Czy Ness

0

Najszybsza metoda, jaką do tej pory wypróbowałem, opiera się na funkcji książki kucharskiej 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

Zobacz odpowiedź, aby uzyskać wyjaśnienie przyspieszenia.


0

Mogę się spóźnić na imprezę, ale będę musiał do tego dodać własny kod. Wykorzystuje w przybliżeniu n / 2 w przestrzeni, ponieważ nie musimy przechowywać liczb parzystych, a także korzystam z modułu python bitarray, dodatkowo drastycznie zmniejszając zużycie pamięci i umożliwiając obliczenie wszystkich liczb pierwszych do 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

Zostało to uruchomione na 64-bitowym 2.4GHZ MAC OSX 10.8.3


1
opublikowanie jednego momentu dla nieznanej maszyny nic nie mówi. Przyjęta odpowiedź brzmi „bez psyco, dla n = 1000000, rwh_primes2 był najszybszy”. Więc jeśli chcesz podać swoje czasy dla tego kodu, jak twoja, na tej samej maszynie, a na 2, 4, 10 mln, jak również, wtedy to byłoby o wiele bardziej pouczające.
Czy Ness

-1, ten kod zależy od specjalnych cech bitrray zaimplementowanych w C, dlatego kod jest szybki, ponieważ większość pracy jest wykonywana w natywnym kodzie w przypisaniu wycinka. Pakiet bitarray WAKACJE standardową definicję właściwego plasterki (indeksowane w pewnym zakresie) do sekwencji modyfikowalnych w to, że umożliwia przypisanie jednego logiczną 0/1 lub Prawda / Fałsz do wszystkich elementów w kawałkach, natomiast średnia dla zachowania czystego Python wydaje się nie zezwalaj na to i zezwalaj tylko na wartość przypisania 0, w którym to przypadku jest to traktowane jako del wszystkich elementów wycinka z sekwencji / tablicy.
GordonBGood,

ciąg dalszy: Jeśli porównano by wywoływanie niestandardowego kodu natywnego, równie dobrze możemy napisać pakiet generatora sekwencji „fastprimes” oparty na kodzie C, takim jak kod pierwszej wersji Kim Walisch, i wygenerować wszystkie liczby pierwsze w czterech miliardach plus 32 -bitowy zakres liczb w ciągu kilku sekund z jednym wywołaniem do generatora sekwencji. Nie wykorzystałoby to prawie żadnej pamięci, ponieważ połączony kod jest oparty na podzielonym na segmenty sicie Eratostenesa, a zatem wykorzystuje tylko kilka dziesiątek kilobajtów pamięci RAM, a gdyby wygenerowano sekwencję, nie byłoby wymagane przechowywanie listy.
GordonBGood,

0

Z czasem zebrałem kilka sit liczb pierwszych. Najszybszy na moim komputerze jest to:

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

Powoli odpowiadam na to pytanie, ale wydawało mi się, że to zabawne ćwiczenie. Używam numpy, który może oszukiwać i wątpię, aby ta metoda była najszybsza, ale powinna być jasna. Przesiewa tablicę boolowską odnoszącą się tylko do jej wskaźników i wywołuje liczby pierwsze ze wskaźników wszystkich wartości True. Nie wymaga modulo.

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]

jest niepoprawny, np. ajs_primes3a(10)-> array([2, 3, 5, 7, 9]). 9nie jest liczbą pierwszą
jfs

Zauważyłeś skrzynkę, której nie widziałem - dobra robota! Problem polegał na tym, że dla idx w zakresie (3, int (do ** ** 0,5), 2): „co powinno być” dla idx w zakresie (3, int (do ** 0,5) + 1, 2): ”. Dzięki, ale teraz działa.
Alan James Salmoni

Powodem było to, że pętla idx wzrosła do „upto ** 05”, co dla przypadków do 15. włącznie włącznie. Od 16 roku działa dobrze. To był zestaw przypadków, dla których nie testowałem. Dodanie 1 oznacza, że ​​powinno działać dla wszystkich liczb.
Alan James Salmoni

Wydaje się, że teraz działa. Jest to najwolniejsze spośród numpyrozwiązań bazujących na zwracaniu tablicy. Uwaga: żadna prawdziwa implementacja Sieve of Eratosthenes nie używa modulo - nie trzeba o tym wspominać. Możesz użyć mat[idx*idx::idx]zamiast mat[idx*2::idx]. I np.nonzero(mat)[0]zamiast np.where(mat == True)[0].
jfs

Dzięki JF. Testowałem przeciwko prime6 () i uzyskałem wynik szybciej do (IIRC) około 250k, kiedy przejęła go prime6 (). primesfrom2to () był szybszy. Na odległości do 20 m, ajs_primes3a () zajął 0,034744977951ms, prime6 () zajął 0,0222899913788ms, a primesfrom2to () zajął 0,0104751586914ms (ta sama maszyna, to samo obciążenie, najlepszy z 10 taktowań). Jest szczerze lepiej, niż się spodziewałem!
Alan James Salmoni

0

Oto interesująca technika generowania liczb pierwszych (jeszcze nie najskuteczniejszych) przy użyciu listowej interpretacji Pythona:

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]

Przykład i kilka wyjaśnień znajdziesz tutaj

Korzystając z naszej strony potwierdzasz, że przeczytałeś(-aś) i rozumiesz nasze zasady używania plików cookie i zasady ochrony prywatności.
Licensed under cc by-sa 3.0 with attribution required.