Schnellster Weg, um alle Primzahlen unter N aufzulisten

357

Dies ist der beste Algorithmus, den ich finden konnte.

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

Kann es noch schneller gemacht werden?

Dieser Code weist einen Fehler auf: Da numberses sich um einen ungeordneten Satz handelt, gibt es keine Garantie dafür, dass numbers.pop()die niedrigste Zahl aus dem Satz entfernt wird. Trotzdem funktioniert es (zumindest für mich) für einige Eingabenummern:

>>> 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
python math optimization primes jbochi
quelle
Das fragliche Code-Sniplet ist viel schneller, wenn Zahlen wie Zahlen = gesetzt deklariert werden (Bereich (n, 2, -2)). Aber kann Sundaram3 nicht schlagen. Danke für die Frage.
Shekhar
3
Es wäre schön, wenn die Antworten Python 3-Versionen der Funktionen enthalten könnten.
Michael Foukarakis
Sicher gibt es eine Bibliothek, um dies zu tun, damit wir nicht unsere eigene> xkcd rollen müssen, versprochen, dass Python so einfach ist wie import antigravity. Gibt es nicht so etwas wie require 'prime'; Prime.take(10)(Ruby)?
Colonel Panic
2
@ColonelPanic Zufällig habe ich github.com/jaredks/pyprimesieve für Py3 aktualisiert und PyPi hinzugefügt. Es ist sicherlich schneller als diese, aber nicht um Größenordnungen - eher ~ 5x schneller als die besten Numpy-Versionen.
Jared
3
@ColonelPanic: Ich denke, dass es angemessen ist, alte Antworten zu bearbeiten, um festzustellen, dass sie gealtert sind, da dies eine nützlichere Ressource darstellt. Wenn die "akzeptierte" Antwort nicht mehr die beste ist, bearbeiten Sie möglicherweise eine Notiz in der Frage mit einem Update von 2015, um die Leute auf die aktuell beste Methode hinzuweisen.
Peter Cordes

Antworten:

366

Warnung: Die timeit Ergebnisse können aufgrund von Hardware- oder Python-Unterschieden variieren.

Unten finden Sie ein Skript, das eine Reihe von Implementierungen vergleicht:

Vielen Dank an stephan für sieve_wheel_30 , um meine Aufmerksamkeit zu bringen. Der Kredit geht an Robert William Hanks für Primzahlen von 2 bis, Primzahlen von 3 bis, rwh_primes, rwh_primes1 und rwh_primes2.

Von den einfachen Python-Methoden, die mit psyco für n = 1000000 getestet wurden , war rwh_primes1 die am schnellsten getestete.

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

Von den einfachen Python-Methoden, die ohne Psyco für n = 1000000 getestet wurden , war rwh_primes2 die schnellste.

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

Von allen getesteten Methoden, die numpy für n = 1000000 zuließen, war primesfrom2to die am schnellsten getestete.

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

Die Timings wurden mit dem folgenden Befehl gemessen:

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

mit {method}durch jeden der Methodennamen ersetzt.

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: <[email protected]>, 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, <[email protected]>, 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)

Durch Ausführen des Skripts wird getestet, dass alle Implementierungen dasselbe Ergebnis liefern.

unutbu
quelle
4
Wenn Sie an nicht reinem Python-Code interessiert sind, sollten Sie dies überprüfen gmpy- er unterstützt Primzahlen über die next_primeMethode seines mpzTyps ziemlich gut .
Alex Martelli
1
Wenn Sie Pypy verwenden, scheinen diese Benchmarks (die Psyco-Benchmarks) ziemlich falsch zu sein. Überraschenderweise fand ich sieveOfEratosthenes und ambi_sieve_plain am schnellsten mit Pypy. Dies ist, was ich für die nicht numpy gefunden habe gist.github.com/5bf466bb1ee9e5726a52
Ehsan Kia
1
Wenn sich jemand fragt, wie sich die Funktionen hier gegen PG7.8 von Wikibooks für reines Python ohne Psyco oder Pypy verhalten: für n = 1000000: PG7.8: 4.93 s pro Schleife; rwh_primes1: 69 ms pro Schleife; rwh_primes2: 57,1 ms pro Schleife
gaborous
8
Können Sie dies mit PyPy aktualisieren, nachdem Psyco tot ist und PyPy es abgelöst hat?
noɥʇʎԀʎzɐɹƆ
3
Wäre toll, wenn diese Funktionen und Timings für Python3 aktualisiert werden könnten.
CS95
135

Schneller und speichertechnischer reiner Python-Code:

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

oder beginnend mit einem halben Sieb

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

Schneller und speichertechnischer Numpy-Code:

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

Eine schnellere Variante, die mit einem Drittel eines Siebs beginnt:

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

Eine (schwer zu codierende) reine Python-Version des obigen Codes wäre:

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

Leider übernimmt Pure-Python nicht die einfachere und schnellere Art der Zuweisung, und das Aufrufen len()innerhalb der Schleife wie in [False]*len(sieve[((k*k)//3)::2*k])ist zu langsam. Also musste ich improvisieren, um Eingaben zu korrigieren (und mehr Mathe zu vermeiden) und extreme (und schmerzhafte) Mathe-Magie zu machen.

Persönlich finde ich es schade, dass numpy (das so weit verbreitet ist) nicht Teil der Python-Standardbibliothek ist und dass die Verbesserungen in Syntax und Geschwindigkeit von Python-Entwicklern völlig übersehen werden.

Robert William Hanks
quelle
2
Numpy ist jetzt mit Python 3 kompatibel. Die Tatsache, dass es nicht in der Standardbibliothek enthalten ist, ist gut, so dass sie ihren eigenen Veröffentlichungszyklus haben können.
Adam
Um nur Binärwerte in einem Array zu speichern, schlage ich vor bitarray- wie hier verwendet (für das einfachste Hauptsieb
Hiro-Protagonist
Sollte primesfrom2to()sich die Teilung beim Gießen in der Klammer innerhalb der Klammern befinden?
355durch113
3
Für eine reine Python-Version, die mit Python 3 kompatibel ist, folgen Sie diesem Link: stackoverflow.com/a/33356284/2482582
Moebius
1
Holy Buttsnacks, dieser Trottel ist schnell.
Scott
42

Es ist eine ziemlich saubere Probe aus dem Python - Kochbuch hier - die schnellste Version auf dieser URL vorgeschlagen ist:

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

das würde also geben

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

Beim Messen an der Shell-Eingabeaufforderung (wie ich es vorziehen würde) mit diesem Code in pri.py beobachte ich:

$ 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

Es sieht also so aus, als ob die Kochbuchlösung doppelt so schnell ist.

Alex Martelli
quelle
1
@jbochi, Sie sind willkommen - aber sehen Sie sich diese URL an, einschließlich der Credits: Zehn von uns haben den Code bis zu diesem Punkt gemeinsam verfeinert, einschließlich Python-Performance-Leuchten wie Tim Peters und Raymond Hettinger (ich habe die geschrieben) endgültiger Text des Rezepts, seit ich das gedruckte Kochbuch bearbeitet habe, aber in Bezug auf die Codierung war mein Beitrag mit dem der anderen vergleichbar ') - am Ende ist es wirklich subtiler und fein abgestimmter Code, und das ist nicht überraschend! -)
Alex Martelli
@Alex: Zu wissen, dass Ihr Code "nur" doppelt so schnell ist wie meiner, macht mich dann ziemlich stolz. :) Die URL war auch sehr interessant zu lesen. Danke noch einmal.
Jbochi
Und es kann mit einer kleinen Änderung noch schneller gemacht werden: siehe stackoverflow.com/questions/2211990/…
tzot
1
... und es kann noch schneller gemacht werden mit einer zusätzlichen ~ 1,2x-1,3x-Beschleunigung, einer drastischen Reduzierung des Speicherbedarfs von O (n) auf O (sqrt (n)) und einer Verbesserung der empirischen Zeitkomplexität, indem die Hinzufügung von verschoben wird Primzahlen für das Diktat, bis ihr Quadrat in der Eingabe zu sehen ist. Testen Sie es hier .
Will Ness
28

Ich glaube, ich habe mit Sundarams Sieb den Rekord von Pure-Python gebrochen:

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)

Vergleich:

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
quelle
1
Ich habe es geschafft, Ihre Funktion um etwa 20% zu beschleunigen, indem ich oben in der Funktion "Null = 0" hinzugefügt und dann das Lambda in Ihrem Filter durch "Null .__ sub__" ersetzt habe. Nicht der schönste Code der Welt, aber etwas schneller :)
Truppo
1
@truppo: Danke für deinen Kommentar! Ich habe gerade festgestellt, dass das Übergeben Noneanstelle der ursprünglichen Funktion funktioniert und sogar noch schneller ist alszero.__sub__
jbochi
7
Wussten Sie, dass wenn Sie bestehen, sundaram3(9)es zurückkehren wird [2, 3, 5, 7, 9]? Es scheint dies mit zahlreichen - vielleicht allen - ungeraden Zahlen zu tun (auch wenn sie keine Primzahlen sind)
wrhall
1
es hat ein Problem: Sundaram3 (7071) enthält 7071, während es nicht Prime ist
bigOther
18

Der Algorithmus ist schnell, hat aber einen schwerwiegenden Fehler:

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

Sie gehen davon aus, dass numbers.pop()dies die kleinste Zahl im Satz zurückgeben würde, dies ist jedoch überhaupt nicht garantiert. Mengen sind ungeordnet und pop()entfernen und geben ein beliebiges Element zurück, sodass es nicht verwendet werden kann, um die nächste Primzahl aus den verbleibenden Zahlen auszuwählen.

etw
quelle
17

Für eine wirklich schnellste Lösung mit ausreichend großem N wäre es, eine vorberechnete Liste von Primzahlen herunterzuladen , als Tupel zu speichern und Folgendes zu tun:

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

Wenn N > primes[-1] nur dann mehr Primzahlen berechnet werden und die neue Liste in Ihrem Code gespeichert wird, ist sie beim nächsten Mal genauso schnell.

Denken Sie immer über den Tellerrand hinaus.

Kimvais
quelle
9
Um fair zu sein, müssten Sie jedoch die Zeit zum Herunterladen, Entpacken und Formatieren der Primzahlen zählen und diese mit der Zeit zum Generieren von Primzahlen mithilfe eines Algorithmus vergleichen. Jeder dieser Algorithmen kann die Ergebnisse problemlos für später in eine Datei schreiben verwenden. Ich denke, in diesem Fall wäre die Numpy-Lösung immer noch schneller, wenn genügend Speicher vorhanden wäre, um tatsächlich alle Primzahlen unter 982.451.653 zu berechnen.
Daniel G
3
@ Daniel richtig. Jedoch steht der Laden, den Sie haben und fahren Sie fort, wann immer nötig, noch ...
Kimvais
@ Daniel GI denke, Download-Zeit ist irrelevant. Geht es nicht wirklich darum, die Zahlen zu generieren? Sie sollten also den Algorithmus berücksichtigen, der zum Erstellen der heruntergeladenen Liste verwendet wird. Und jede Zeitkomplexität würde die einmalige Dateiübertragung bei O (n) ignorieren.
Ross
Die FAQ für die UTM - Prime Seite schlägt vor , die Berechnung kleine Primzahlen ist schneller als sie aus einer Platte zu lesen (die Frage ist , was kleine Mittel).
Batman
12

Wenn Sie das Rad nicht neu erfinden möchten, können Sie die symbolische Mathematikbibliothek sympy installieren (ja, es ist Python 3-kompatibel).

pip install sympy

Und verwenden Sie die Primerange- Funktion

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

Wenn Sie itertools akzeptieren, aber nicht numpy, finden Sie hier eine Anpassung von rwh_primes2 für Python 3, die auf meinem Computer etwa doppelt so schnell ausgeführt wird. Die einzige wesentliche Änderung besteht darin, ein Bytearray anstelle einer Liste für den Booleschen Wert zu verwenden und die Komprimierung anstelle eines Listenverständnisses zu verwenden, um die endgültige Liste zu erstellen. (Ich würde dies als Kommentar wie moarningsun hinzufügen, wenn ich könnte.)

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

Vergleiche:

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

und

>>> 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
Jason
quelle
Sehr coole Implementierung. :)
Krish
7

Es ist lehrreich, Ihren eigenen Prime-Finding-Code zu schreiben, aber es ist auch nützlich, eine schnelle, zuverlässige Bibliothek zur Hand zu haben. Ich habe einen Wrapper um das Primesieve der C ++ - Bibliothek geschrieben , der Primesieve-Python heißt

Versuch es pip install primesieve

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

Ich wäre gespannt auf die Geschwindigkeit im Vergleich.

Oberst Panik
quelle
Es ist nicht genau das, was OP bestellt hat, aber ich verstehe nicht, warum die Abstimmung. Es ist eine 2,8-Sekunden-Lösung im Gegensatz zu einigen anderen externen Modulen. Ich habe in der Quelle bemerkt, dass es Threaded ist, habe irgendwelche Tests, wie gut es skaliert?
Jetibo
@ Jetibo Prost. Der Engpass scheint darin zu bestehen, den C ++ - Vektor in die Python-Liste zu kopieren, daher ist die count_primesFunktion viel schneller alsgenerate_primes
Colonel Panic
Auf meinem Computer kann es bequem Primzahlen bis zu 1e8 erzeugen (es gibt MemoryError für 1e9) und Primzahlen bis zu 1e10 zählen. @HappyLeapSecond oben vergleicht Algorithmen für 1e6
Colonel Panic
7

Hier sind zwei aktualisierte (reine Python 3.6) Versionen einer der schnellsten Funktionen:

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:])]
Bruno Astrolino
quelle
1
In Python 3 habe ich diese Funktion stackoverflow.com/a/3035188/7799269 verwendet, aber / durch // und xrange durch range ersetzt und sie schienen viel schneller zu sein.
samerivertwice
4

Eine deterministische Implementierung des Miller-Rabin-Primalitätstests unter der Annahme, dass N <9.080.191 ist

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)

Laut dem Artikel auf Wikipedia ( http://en.wikipedia.org/wiki/Miller–Rabin_primality_test ) reicht es aus, N <9.080.191 auf a = 2,3,37 zu testen, und 73 reicht aus, um zu entscheiden, ob N zusammengesetzt ist oder nicht.

Und ich habe den Quellcode aus der probabilistischen Implementierung des ursprünglichen Miller-Rabin-Tests angepasst, der hier zu finden ist: http://en.literateprograms.org/Miller-Rabin_primality_test_(Python)

Ruggiero Spearman
quelle
1
Vielen Dank für den Miller-Rabin-Primalitätstest, aber dieser Code ist tatsächlich langsamer und liefert nicht die richtigen Ergebnisse. 37 ist prim und besteht den Test nicht.
Jbochi
Ich denke, 37 ist einer der Sonderfälle, mein schlechtes. Ich war hoffnungsvoll über die deterministische Version :)
Ruggiero Spearman
Es gibt keinen Sonderfall für Rabin Miller.
Fehlgeleitet
2
Sie haben den Artikel falsch gelesen. Es ist 31, nicht 37. Aus diesem Grund schlägt Ihre Implementierung fehl.
Logan
4

Wenn Sie die Kontrolle über N haben, können Sie alle Primzahlen am schnellsten auflisten, indem Sie sie vorberechnen. Ernsthaft. Precomputing ist eine Möglichkeit, die Optimierung zu übersehen.

Dave W. Smith
quelle
3
Oder laden Sie sie von hier herunter. Primes.utm.edu/lists/small/millions , aber die Idee ist, das Python-Limit zu testen und zu sehen, ob aus der Optimierung schöner Code hervorgeht.
Jbochi
4

Hier ist der Code, den ich normalerweise zum Generieren von Primzahlen in Python verwende:

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

Es kann nicht mit den hier veröffentlichten schnelleren Lösungen mithalten, aber es ist zumindest reine Python.

Vielen Dank für die Veröffentlichung dieser Frage. Ich habe heute wirklich viel gelernt.

MAK
quelle
3

Für den schnellsten Code ist die Numpy-Lösung die beste. Aus rein akademischen Gründen veröffentliche ich jedoch meine reine Python-Version, die etwas weniger als 50% schneller ist als die oben veröffentlichte Kochbuchversion. Da ich die gesamte Liste im Speicher erstelle, benötigen Sie genügend Speicherplatz, um alles aufzunehmen, aber es scheint ziemlich gut zu skalieren.

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)

Und die Ergebnisse:

>>>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
Daniel G.
quelle
3

Eine etwas andere Implementierung eines halben Siebs mit Numpy:

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

Mathe importieren
numpy importieren
def prime6 (bis):
    Primzahlen = numpy.arange (3, bis zu + 1,2)
    isprime = numpy.ones ((bis zu 1) / 2, dtype = bool)
    für Faktor in Primzahlen [: int (math.sqrt (bis zu))]:
        wenn isprime [(Faktor-2) / 2]: isprime [(Faktor * 3-2) / 2: (bis zu 1) / 2: Faktor] = 0
    return numpy.insert (Primzahlen [isprime], 0,2)

Kann jemand dies mit den anderen Timings vergleichen? Auf meiner Maschine scheint es ziemlich vergleichbar mit dem anderen Numpy-Halbsieb zu sein.

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

Es ist alles geschrieben und getestet. Das Rad muss also nicht neu erfunden werden.

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

gibt uns einen Rekord von 12,2 ms !

10 loops, best of 10: 12.2 msec per loop

Wenn dies nicht schnell genug ist, können Sie PyPy ausprobieren:

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

was in ... resultiert:

10 loops, best of 10: 2.03 msec per loop

Die Antwort mit 247 Up-Votes listet 15,9 ms für die beste Lösung auf. Vergleichen Sie das !!!

lifolofi
quelle
3

Ich habe einige Funktionen von unutbu getestet und sie mit hungrigen Millionen berechnet

Die Gewinner sind die Funktionen, die die Numpy-Bibliothek verwenden.

Hinweis : Es wäre auch interessant, einen Speicherauslastungstest durchzuführen :)

Ergebnis der Rechenzeit

Beispielcode

Vollständiger Code in meinem Github-Repository

#!/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()
Bruno Adelé
quelle
2
Um algorithmische Leistungen zu vergleichen , ist es besser, auf einer Log-Log-Skala zu zeichnen .
Will Ness
3

Für 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
quelle
3

Schnellstes Hauptsieb in 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

Ich habe Sieb of Eratosthenes für Geschwindigkeit und Gedächtnis optimiert .

Benchmark

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)

Ausgabe

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

Wenn Sie Python zum ersten Mal verwenden, scheinen einige der Methoden, die ich hier verwende, etwas umständlich zu sein. Ich habe gerade meinen C ++ - Code direkt in Python konvertiert und das ist es, was ich habe (wenn auch ein bisschen langsam in 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

Gefunden 664579 Primzahlen in 12.799119 Sekunden!

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

Gefunden 664579 Primzahlen in 10.230172 Sekunden!

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

Gefunden 664579 Primzahlen in 7.113776 Sekunden!

smac89
quelle
2

Ich weiß, dass der Wettbewerb seit einigen Jahren geschlossen ist. …

Dies ist jedoch mein Vorschlag für ein reines Python-Hauptsieb, bei dem die Vielfachen von 2, 3 und 5 weggelassen werden, indem geeignete Schritte ausgeführt werden, während das Sieb vorwärts verarbeitet wird. Trotzdem ist es für N <10 ^ 9 tatsächlich langsamer als für @Robert William Hanks überlegene Lösungen rwh_primes2 und rwh_primes1. Durch die Verwendung eines ctypes.c_ushort-Siebarrays über 1,5 * 10 ^ 8 ist es irgendwie an Speichergrenzen anpassbar.

10 ^ 6

$ python -mtimeit -s "importiere primeSieveSpeedComp" "primeSieveSpeedComp.primeSieveSeq (1000000)" 10 Schleifen, am besten 3: 46,7 ms pro Schleife

zum Vergleichen: $ python -mtimeit -s "importiere primeSieveSpeedComp" "primeSieveSpeedComp.rwh_primes1 (1000000)" 10 Schleifen, am besten 3: 43,2 ms pro Schleife zum Vergleichen: $ python -m timeit -s "importiere primeSieveSpeedComp" "primeSieveSpeedComp.rwh_ (1000000) 10 Schleifen, am besten 3: 34,5 ms pro Schleife

10 ^ 7

$ python -mtimeit -s "importiere primeSieveSpeedComp" "primeSieveSpeedComp.primeSieveSeq (10000000)" 10 Schleifen, am besten 3: 530 ms pro Schleife

zum Vergleichen: $ python -mtimeit -s "importiere primeSieveSpeedComp" "primeSieveSpeedComp.rwh_primes1 (10000000)" 10 Schleifen, am besten 3: 494 ms pro Schleife zum Vergleichen: $ python -m timeit -s "importiere primeSieveSpeedComp" "primeSieveSpeedComp.rwh_ (10000000) 10 Schleifen, am besten 3: 375 ms pro Schleife

10 ^ 8

$ python -mtimeit -s "importiere primeSieveSpeedComp" "primeSieveSpeedComp.primeSieveSeq (100000000)" 10 Schleifen, am besten 3: 5,55 Sek. pro Schleife

zum Vergleichen: $ python -mtimeit -s "import primeSieveSpeedComp" "primeSieveSpeedComp.rwh_primes1 (100000000)" 10 Schleifen, am besten 3: 5,33 Sekunden pro Schleife zum Vergleichen: $ python -m timeit -s "import primeSieveSpeedComp" "primeSieveSpeedComp.rwh_ (100000000) 10 Schleifen, am besten 3: 3,95 Sekunden pro Schleife

10 ^ 9

$ python -mtimeit -s "importiere primeSieveSpeedComp" "primeSieveSpeedComp.primeSieveSeq (1000000000)" 10 Schleifen, am besten 3: 61,2 Sekunden pro Schleife

Zum Vergleich: $ python -mtimeit -n 3 -s "import primeSieveSpeedComp" "primeSieveSpeedComp.rwh_primes1 (1000000000)" 3 Schleifen, am besten 3: 97,8 Sekunden pro Schleife

Zum Vergleich: $ python -m timeit -s "importiere primeSieveSpeedComp" "primeSieveSpeedComp.rwh_primes2 (1000000000)" 10 Schleifen, am besten 3: 41,9 Sekunden pro Schleife

Sie können den folgenden Code in ubuntus primeSieveSpeedComp kopieren, um diese Tests zu überprüfen.

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
quelle
Um Ihre Testergebnisse zu visualisieren, zeichnen Sie sie auf einer Log-Log-Skala auf, um die empirischen Wachstumsordnungen zu sehen und zu vergleichen .
Will Ness
@ Will danke für die Eingabe, ich werde dies berücksichtigen, wenn ich das nächste Mal einen solchen Vergleich brauche
ABri
1

Hier ist eine numpy-Version von Sieve of Eratosthenes, die sowohl eine gute Komplexität (niedriger als das Sortieren eines Arrays der Länge n) als auch eine Vektorisierung aufweist. Im Vergleich zu @unutbu mal ist dies genauso schnell wie die Pakete mit 46 Mikrosekunden, um alle Primzahlen unter einer Million zu finden.

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]

Timings: 

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
Peter Mølgaard Pallesen
quelle
1

Ich habe einen Großteil des Codes für Python 3 aktualisiert und ihn auf perfplot (ein Projekt von mir) geworfen, um zu sehen, welcher tatsächlich am schnellsten ist. Es stellt sich heraus, dass für große n, primesfrom{2,3}tonehmen Sie den Kuchen:

Geben Sie hier die Bildbeschreibung ein

Code zur Reproduktion der Handlung:

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: <[email protected]>, 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, <[email protected]>, 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",
)
Nico Schlömer
quelle
0

Ich vermute, dass der schnellste Weg darin besteht, die Primzahlen in Ihrem Code hart zu codieren.

Warum also nicht einfach ein langsames Skript schreiben, das eine andere Quelldatei generiert, in der alle Zahlen fest verdrahtet sind, und diese Quelldatei dann importieren, wenn Sie Ihr eigentliches Programm ausführen?

Dies funktioniert natürlich nur, wenn Sie die Obergrenze von N zur Kompilierungszeit kennen, ist jedoch bei (fast) allen Projekt-Euler-Problemen der Fall.

 

PS: Ich könnte mich zwar irren, wenn das Parsen der Quelle mit fest verdrahteten Primzahlen langsamer ist als das Berechnen, aber soweit ich weiß, läuft Python aus kompilierten .pycDateien, daher sollte das Lesen eines binären Arrays mit allen Primzahlen bis zu N blutig sein schnell in diesem Fall.

akuhn
quelle
0

Tut mir leid, aber erat2 () hat einen schwerwiegenden Fehler im Algorithmus.

Bei der Suche nach dem nächsten Verbund müssen wir nur ungerade Zahlen testen. q, p beide sind ungerade; dann ist q + p gerade und muss nicht getestet werden, aber q + 2 * p ist immer ungerade. Dies eliminiert den "wenn gerade" -Test in der while-Schleifenbedingung und spart etwa 30% der Laufzeit.

Während wir gerade dabei sind: Verwenden Sie anstelle der eleganten Methode 'D.pop (q, None)' get and delete 'if q in D: p = D [q], del D [q]', was doppelt so schnell ist ! Zumindest auf meiner Maschine (P3-1Ghz). Daher schlage ich diese Implementierung dieses cleveren Algorithmus vor:

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
quelle
Eine verschobene Hinzufügung von Primzahlen zum Diktat (bis das Quadrat einer Primzahl in der Eingabe angezeigt wird ) finden Sie unter stackoverflow.com/a/10733621/849891 .
Will Ness
0

Die schnellste Methode, die ich bisher ausprobiert habe, basiert auf der Python-Kochbuchfunktionerat2 :

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

In dieser Antwort finden Sie eine Erklärung für die Beschleunigung.

tzot
quelle
0

Ich komme möglicherweise zu spät zur Party, muss aber dafür meinen eigenen Code hinzufügen. Es benötigt ungefähr n / 2 im Speicherplatz, da wir keine geraden Zahlen speichern müssen, und ich verwende auch das Bitarray-Python-Modul, um den Speicherverbrauch weiter drastisch zu reduzieren und die Berechnung aller Primzahlen bis zu 1.000.000.000 zu ermöglichen

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

Dies wurde auf einem 64-Bit-2,4-GHz-MAC OSX 10.8.3 ausgeführt

Cobie
quelle
1
Das Posten eines Timings für eine unbekannte Maschine sagt nichts aus. Die hier akzeptierte Antwort lautet "ohne Psyco war rwh_primes2 für n = 1000000 am schnellsten". Also , wenn Sie Ihre Timings für diesen Code als auch bei Ihnen zur Verfügung stellen würde, auf der gleichen Maschine, und bei 2, 4, 10 Millionen als gut, dann würde es viel mehr informativ.
Will Ness
-1, Dieser Code hängt von den speziellen Funktionen des in C implementierten Bitarrays ab, weshalb der Code schnell ist, da der größte Teil der Arbeit in nativem Code in der Slice-Zuweisung ausgeführt wird. Das Bitarray-Paket bricht die Standarddefinition für richtige Slices (über einen Bereich indiziert) für veränderbare Sequenzen, indem es die Zuweisung eines einzelnen Booleschen Wertes 0/1 oder True / False zu allen Elementen des Slice ermöglicht, während das Standardverhalten für reines Python zu sein scheint Dies nicht zulassen und nur den Zuweisungswert 0 zulassen. In diesem Fall wird er als del aller Slice-Elemente aus der Sequenz / dem Array behandelt.
GordonBGood
Fortsetzung: Wenn der Aufruf von nicht standardmäßigem nativem Code verglichen werden soll, können wir auch ein "Fastprimes" -Sequenzgeneratorpaket schreiben, das auf C-Code wie dem von Kim Walischs Primesieve basiert, und alle Primzahlen in den vier Milliarden plus 32 generieren -Bit-Nummernbereich in nur wenigen Sekunden mit einem einzigen Aufruf des Sequenzgenerators. Dies würde auch fast keinen Speicher verbrauchen, da der verknüpfte Code auf einem segmentierten Sieb von Eratosthenes basiert und daher nur einige Zehntel Kilobyte RAM verwendet, und wenn eine Sequenz erzeugt würde, wäre kein Listenspeicher erforderlich.
GordonBGood
0

Ich habe im Laufe der Zeit mehrere Primzahlsiebe gesammelt. Das schnellste auf meinem Computer ist das Folgende:

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
Stefan Gruenwald
quelle
0

Ich beantworte diese Frage nur langsam, aber es schien eine lustige Übung zu sein. Ich benutze Numpy, was betrügen könnte und ich bezweifle, dass diese Methode die schnellste ist, aber es sollte klar sein. Es siebt ein Boolesches Array, das sich nur auf seine Indizes bezieht, und ermittelt Primzahlen aus den Indizes aller True-Werte. Kein Modulo erforderlich.

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]
Alan James Salmoni
quelle
es ist falsch zB ajs_primes3a(10)-> array([2, 3, 5, 7, 9]).9ist kein Prime
JFS
Sie haben einen Randfall entdeckt, den ich nicht hatte - gut gemacht! Das Problem war in 'für idx im Bereich (3, int (bis zu ** 0,5), 2):' was sein sollte 'für idx im Bereich (3, int (bis zu ** 0,5) + 1, 2):'. Danke, aber es funktioniert jetzt.
Alan James Salmoni
Der Grund war, dass die idx-Schleife auf 'bis ** 05' stieg, was für Fälle bis einschließlich 15 gilt. Ab 16 funktioniert sie einwandfrei. Dies war eine Reihe von Randfällen, auf die ich nicht getestet hatte. Das Hinzufügen von 1 bedeutet, dass es für alle Zahlen funktionieren sollte.
Alan James Salmoni
Es scheint jetzt zu funktionieren. Es ist die langsamste unter den numpybasierten Lösungen, die ein Array zurückgeben. Hinweis: Keine echte Sieve of Eratosthenes-Implementierung verwendet Modulo - es muss nicht erwähnt werden. Sie könnten mat[idx*idx::idx]anstelle von verwenden mat[idx*2::idx]. Und np.nonzero(mat)[0]statt np.where(mat == True)[0].
JFS
Danke JF. Ich habe gegen prime6 () getestet und ein schnelleres Ergebnis bis zu (IIRC) von ca. 250k erzielt, als prime6 () übernahm. primesfrom2to () war schneller. Bei bis zu 20 m dauerte ajs_primes3a () 0,034744977951 ms, prime6 () 0,0222899913788 ms und primesfrom2to () 0,0104751586914 ms (gleiche Maschine, gleiche Last, beste von 10 Timings). Es ist ehrlich gesagt besser als ich dachte!
Alan James Salmoni
0

Hier ist eine interessante Technik zum Generieren von Primzahlen (aber nicht die effizienteste) unter Verwendung des Listenverständnisses von 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]

Das Beispiel und einige Erklärungen finden Sie hier

Alexander
quelle