Schnellster Weg, um alle Primzahlen unter N aufzulisten

Question

Dies ist der beste Algorithmus, den ich entwickeln 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 er noch schneller gemacht werden?

Dieser Code hat ein Problem: Da numbers eine ungeordnete Menge ist, gibt es keine Garantie dafür, dass numbers.pop() die niedrigste Zahl aus der Menge entfernt. Trotzdem funktioniert es (zumindest für mich) für einige Eingabelnummern:

>>> sum(get_primes(2000000))
142913828922L
#Das ist die korrekte Summe aller Zahlen unter 2 Millionen
>>> 529 in get_primes(1000)
False
>>> 529 in get_primes(530)
True

Answer 1

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

Dies ist dennoch mein Vorschlag für ein reines Python-Primzahl-Sieb, basierend auf der Auslassung der Vielfachen von 2, 3 und 5 durch die Verwendung geeigneter Schritte beim Voranschreiten des Siebes. Trotzdem ist es tatsächlich langsamer für N<10^9 als die überlegenen Lösungen von @Robert William Hanks rwh_primes2 und rwh_primes1. Durch die Verwendung eines ctypes.c_ushort-Siebarrays über 1,5 * 10^8 ist es irgendwie an die Speichergrenzen anpassbar.

10^6

$ python -mtimeit -s"import primeSieveSpeedComp" "primeSieveSpeedComp.primeSieveSeq(1000000)" 10 Schleifen, bester von 3: 46,7 msec pro Schleife

Zum Vergleich: $ python -mtimeit -s"import primeSieveSpeedComp" "primeSieveSpeedComp.rwh_primes1(1000000)" 10 Schleifen, bester von 3: 43,2 msec pro Schleife Zum Vergleich: $ python -m timeit -s"import primeSieveSpeedComp" "primeSieveSpeedComp.rwh_primes2(1000000)" 10 Schleifen, bester von 3: 34,5 msec pro Schleife

10^7

$ python -mtimeit -s"import primeSieveSpeedComp" "primeSieveSpeedComp.primeSieveSeq(10000000)" 10 Schleifen, bester von 3: 530 msec pro Schleife

Zum Vergleich: $ python -mtimeit -s"import primeSieveSpeedComp" "primeSieveSpeedComp.rwh_primes1(10000000)" 10 Schleifen, bester von 3: 494 msec pro Schleife Zum Vergleich: $ python -m timeit -s"import primeSieveSpeedComp" "primeSieveSpeedComp.rwh_primes2(10000000)" 10 Schleifen, bester von 3: 375 msec pro Schleife

10^8

$ python -mtimeit -s"import primeSieveSpeedComp" "primeSieveSpeedComp.primeSieveSeq(100000000)" 10 Schleifen, bester von 3: 5,55 Sek. pro Schleife

Zum Vergleich: $ python -mtimeit -s"import primeSieveSpeedComp" "primeSieveSpeedComp.rwh_primes1(100000000)" 10 Schleifen, bester von 3: 5,33 Sek. pro Schleife Zum Vergleich: $ python -m timeit -s"import primeSieveSpeedComp" "primeSieveSpeedComp.rwh_primes2(100000000)" 10 Schleifen, bester von 3: 3,95 Sek. pro Schleife

10^9

$ python -mtimeit -s"import primeSieveSpeedComp" "primeSieveSpeedComp.primeSieveSeq(1000000000)" 10 Schleifen, bester von 3: 61,2 Sek. pro Schleife

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

Zum Vergleich: $ python -m timeit -s"import primeSieveSpeedComp" "primeSieveSpeedComp.rwh_primes2(1000000000)" 10 Schleifen, bester von 3: 41,9 Sek. pro Schleife

Sie können den untenstehenden Code in dem PrimeSieveSpeedComp in Ubuntu 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 'benutzt ctypes "unsigned short int Array"'
    else:
        sieve = (MAX_Int >> 1) * [False]
        #print 'benutzt python list() von 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

Answer 2

Vielleicht bin ich langsam auf die Party gekommen, aber ich mache es mit dem schnellsten Code wieder wett (zumindest auf meinem Computer, zum Zeitpunkt des Schreibens). Dieser Ansatz verwendet sowohl numpy als auch bitarray und ist inspiriert von primesfrom2to aus dieser Antwort.

import numpy as np
from bitarray import bitarray

def bit_primes(n):
    bit_sieve = bitarray(n // 3 + (n % 6 == 2))
    bit_sieve.setall(1)
    bit_sieve[0] = False

    for i in range(int(n ** 0.5) // 3 + 1):
        if bit_sieve[i]:
            k = 3 * i + 1 | 1
            bit_sieve[k * k // 3::2 * k] = False
            bit_sieve[(k * k + 4 * k - 2 * k * (i & 1)) // 3::2 * k] = False

    np_sieve = np.unpackbits(np.frombuffer(bit_sieve.tobytes(), dtype=np.uint8)).view(bool)
    return np.concatenate(((2, 3), ((3 * np.flatnonzero(np_sieve) + 1) | 1)))

Hier ist ein Vergleich mit primesfrom2to, das zuvor als schnellste Lösung in unutbu's Vergleich gefunden wurde:

python3 -m timeit -s "import fast_primes" "fast_primes.bit_primes(1000000)"
200 loops, best of 5: 1.19 msec per loop

python3 -m timeit -s "import fast_primes" "fast_primes.primesfrom2to(1000000)"
200 loops, best of 5: 1.23 msec per loop

Für das Auffinden von Primzahlen unter 1 Million war bit_primes etwas schneller. Bei größeren Werten von n kann der Unterschied deutlicher sein. In einigen Fällen war bit_primes mehr als doppelt so schnell:

python3 -m timeit -s "import fast_primes" "fast_primes.bit_primes(500_000_000)"
1 loop, best of 5: 540 msec per loop

python3 -m timeit -s "import fast_primes" "fast_primes.primesfrom2to(500_000_000)"
1 loop, best of 5: 1.15 sec per loop

Zu Referenzzwecken hier die minimal modifizierte (für Python 3 angepasste) Version von primesfrom2to, mit der ich verglichen habe:

def primesfrom2to(n):
    # https://stackoverflow.com/questions/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 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)]

Answer 3

Es überrascht mich, dass noch niemand numba erwähnt hat.

Diese Version erreicht die 1-Millionen-Marke in 2,47 ms ± 36,5 µs.

Vor Jahren wurde auf der Wikipedia-Seite Primzahl Pseudocode für eine Version des Atkin-Siebs gegeben. Dies ist dort nicht mehr verfügbar, und eine Referenz zum Atkin-Sieb scheint einen anderen Algorithmus darzustellen. Eine Version der Wikipedia-Seite vom 01.03.2007, Primzahl vom 01.03.2007, zeigt den von mir verwendeten Pseudocode.

import numpy as np
from numba import njit

@njit
def nb_primes(n):
    # Generiert Primzahlen 2 <= p <= n
    # Atkin-Sieb -- siehe https://en.wikipedia.org/w/index.php?title=Prime_number&oldid=111775466
    sqrt_n = int(np.sqrt(n)) + 1

    # initialisiere das Sieb
    s = np.full(n + 1, -1, dtype=np.int8)
    s[2] = 1
    s[3] = 1

    # füge mögliche Primzahlen ein:
    # ganze Zahlen, die eine ungerade Anzahl von
    # Darstellungen durch bestimmte quadratische Formen haben
    for x in range(1, sqrt_n):
        x2 = x * x
        for y in range(1, sqrt_n):
            y2 = y * y
            k = 4 * x2 + y2
            if k <= n and (k % 12 == 1 or k % 12 == 5): s[k] *= -1
            k = 3 * x2 + y2
            if k <= n and (k % 12 == 7): s[k] *= -1
            k = 3 * x2 - y2
            if k <= n and x > y and k % 12 == 11: s[k] *= -1

    # eliminiere zusammengesetzte Zahlen durch Sieben
    for k in range(5, sqrt_n):
        if s[k]:
            k2 = k*k
            # k ist eine Primzahl, lasse Vielfache ihres Quadrats aus; dies genügt, da
            # zusammengesetzte Zahlen, die es auf die Liste geschafft haben, nicht quadratfrei sein können
            for i in range(1, n // k2 + 1):
                j = i * k2 # j {k², 2k², 3k², ..., n}
                s[j] = -1
    return np.nonzero(s>0)[0]

# erste Ausführung für "Kompilierung"
nb_primes(10)

Zeitmessung

In[10]:
%timeit nb_primes(1_000_000)

Out[10]:
2,47 ms ± 36,5 µs pro Durchlauf (Mittelwert ± Standardabweichung von 7 Durchläufen, 100 Durchläufe pro Durchlauf)

In[11]:
%timeit nb_primes(10_000_000)

Out[11]:
33,4 ms ± 373 µs pro Durchlauf (Mittelwert ± Standardabweichung von 7 Durchläufen, 10 Durchläufe pro Durchlauf)

In[12]:
%timeit nb_primes(100_000_000)

Out[12]:
828 ms ± 5,64 ms pro Durchlauf (Mittelwert ± Standardabweichung von 7 Durchläufen, 1 Durchlauf pro Durchlauf)

Answer 4

Für den schnellsten Code ist die numpy-Lösung am besten. Aus rein akademischen Gründen poste ich jedoch meine reine Python-Version, die etwas weniger als 50% schneller ist als die oben gepostete Kochbuch-Version. Da ich die gesamte Liste im Speicher erstelle, benötigen Sie genügend Platz, um alles zu halten, aber es scheint ziemlich gut zu skalieren.

def daniel_sieve_2(maxNumber):
    """
    Bei einer Zahl liefert alle Zahlen zurück, die kleiner oder gleich dieser Zahl und prim sind.
    """
    allNumbers = range(3, maxNumber+1, 2)
    for mIndex, number in enumerate(xrange(3, maxNumber+1, 2)):
        if allNumbers[mIndex] == 0:
            continue
        # setze jetzt alle Vielfachen auf 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 "Meins: {0:0.4f} ms".format(min(mine.repeat(3, 1))*1000)
Meins: 428.9446 ms
>>>print "Vorher bestes {0:0.4f} ms".format(min(prev.repeat(3, 1))*1000)
Vorher bestes 621.3581 ms

Answer 5

Erstmalige Verwendung von Python, daher mögen einige der von mir verwendeten Methoden etwas umständlich erscheinen. Ich habe einfach meinen C++-Code direkt in Python konvertiert und das ist das Ergebnis (obwohl 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 #Dies ist für 2

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

elapsed = (time.clock() - start)
print "\nGefunden", count, "Primzahlen in", elapsed, "Sekunden!\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 #Dies ist für 2

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

elapsed = (time.clock() - start)
print "\nGefunden", count, "Primzahlen in", elapsed, "Sekunden!\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 #Dies ist für 2

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

elapsed = (time.clock() - start)
print "\nGefunden", count, "Primzahlen in", elapsed, "Sekunden!\n"

python Primes2.py

Gefunden 664579 Primzahlen in 7.113776 Sekunden!