2016-05-05 6 views
5

In Ordnung, so bekomme ich das Prinzip hinter dem Sieb von Eratosthenes. Ich habe versucht und größtenteils gescheitert, selbst einen zu schreiben (ich schrieb einen funktionalen Primzahlgenerator, es war einfach nicht effizient oder ein Sieb). Ich denke, ich habe mehr Probleme mit der Mathematik, aber die Programmierung hat mich auch verwirrt.Könnte jemand mir helfen, dieses Sieb von Eratosthenes Skript zu verstehen? Die letzten paar Zeilen haben mich gestoßen

import numpy as np 

def primesfrom2to(n): 
    sieve = np.ones(n/3 + (n%6==2), dtype=np.bool) 
    # print(sieve) for n = 10 returns [True, True, True] 

Ok, schreiben Sie die Fledermaus Ich bin ein wenig verwirrt. Es erzeugt eine Liste von Wahrheitswerten, die fortschreitend als falsch gekennzeichnet werden, da sie als zusammengesetzt gelten. Ich denke, es wird gesagt, dass nicht mehr als 1/3 der Werte, die kleiner als n sind, prim sind, aber was bringt das Hinzufügen einer Modolus-Operation?

sieve[0] = False 

Marken 1 als falsch?

for i in range(int(n**0.5)//3+1): 
     # print(i) for n = 10 returns 0 and 1 on separate lines 

Hier wird der zu prüfende Bereich eingestellt. Keine Zahl hat ein Produkt größer als seine Quadratwurzel. Was ist mit der Teilung von 3+1?

 if sieve[i]: 
      k=3*i+1|1 
      #print(k) for n = 10 returns '5' doesn't change till it adds 7 at n = 50 

Ok so, wenn sieve[i] wahr ist (so Ansaug-/noch nicht als Verbund markiert?), Dann 3*i + 1 or 1? Wie funktioniert das genau? Es scheint Primzahlen zu erzeugen, dass es später multipliziert wird, um alle nachfolgenden Produkte zu entfernen, aber ich weiß nicht wie.

  sieve[  ((k*k)/3)  ::2*k] = False 
      sieve[(k*k+4*k-2*k*(i&1))/3::2*k] = False 

ok so diese beiden nehmen die Primzahlen k und Kennzeichnung aller weiteren Vielfachen von ihnen? wenn n=5 würde das nicht sieve[8.33::10] = False machen? Und der andere wie sieve[41.3::10]? Ich verstehe es einfach nicht.

print(np.r_[2,3,((3*np.nonzero(sieve)[0]+1)|1)]) 

Alles in Ordnung, und schließlich tatsächlich die Liste der Primzahlen zu generieren. Nochmal, was ist los mit der Multiplikation mit drei? Offensichtlich gibt es etwas Grundlegendes zu 3 in diesem Code, das ich einfach nicht verstehe. Auch ich kann das |1 Konzept nicht wieder erfassen.

Oh und nur zum Spaß hier ist meine effektive und schrecklich ineffiziente Versuch.

import numpy 
def sieve(num): 

    master = numpy.array([i for i in range(2, num+1)]) 
    run = [] 
    y=2 


    while y <= ((num+1)**(1/2)): 
     thing = [x for x in range(y, num+1, y) if x > 5 or x == 4] 
     run = run + thing 
     y = y+1 
    alist = numpy.in1d(master, run, invert = True) 
    blist = (master[alist]) 
    print(blist) 

Es dauerte 57s, um die Primzahlen bis zu 500.000 zu berechnen. Ich machte die Eulersumme von Primzahlen bis zu 2.000.000 Problem.

Antwort

3

Hier ist ein einfacher Siebalgorithmus in numpy:

import numpy as np 
def sieve(Nmax): 
    """Calculate prime numbers between 0 and Nmax-1.""" 
    is_prime = np.ones(Nmax, dtype=bool) 
    Pmax = int(np.sqrt(Nmax)+1) 
    is_prime[0] = is_prime[1] = False 
    p = 2 
    dp = 1 
    while p <= Pmax: 
     if is_prime[p]: 
      is_prime[p*2::p] = False 
     p += dp 
     dp = 2 
    primes = np.argwhere(is_prime)[:,0] 
    print(primes) 

sieve(100)  

Wenn ich die print-Anweisung zu entfernen, ist es für Nmax läuft = 10^8 in 2,5 Sekunden auf meinem bescheidenen Laptop.

Dieser Algorithmus ist ein bisschen ineffizient, da er die Primzahl von geraden Werten und Vielfachen von 3 speichern muss. Sie können Speicherplatz sparen, indem Sie diese ausfiltern, sodass Sie die Primzahl von n verfolgen * 6 + 1 und n * 6 + 5, für jede ganze Zahl n> = 1. Das spart Ihnen zwei Drittel des Speicherplatzes, auf Kosten von etwas mehr Buchhaltung. Es scheint, dass Sie versucht haben, mit der schwierigen Version zu beginnen. Außerdem wird die gesamte Buchhaltung teuer, wenn der Python-Interpreter if Anweisungen auswertet und die Speicherverwaltung Ihrer Listen übernimmt. Das Array-Slicing von Python/numpy ist eine viel natürlichere Methode, und es ist wahrscheinlich schnell genug für Ihre Zwecke.

In Bezug auf Ihre Fragen:

  • (n% 6 == 2) boolean, wird als 0 oder 1 interpretiert werden und wird ein weiteres Element füge eine "off-by-one" Fehler zu verhindern.
  • int(n**0.5)//3+1 sollte als int(int(n**0.5)/3) + 1 gelesen werden. Division geht vor Zusatz. Die Division durch drei ist so, dass Sie nur Platz für 6n + 1 und 6n + 5 Werte zuweisen.
  • k=3*i+1|1 bedeutet k=3*i+1, fügen Sie eine hinzu, wenn es gerade ist (es ist ein bisschen weise oder). Das Array-Element 'i' entspricht der potentiellen Primzahl (3*i+1)|1. Wenn also die Array-Indizes [0, 1, 2, 3, 4, 5, ...] sind, repräsentieren sie die Nummern [1, 5, 7, 11, 13, 17, ...].
  • Die Einstellung der Siebelemente auf False erfolgt getrennt für die Elemente 6n + 1 und 6n + 5, deshalb gibt es zwei Zuweisungen.
  • Wenn Sie dies in Python 2.7 ausführen, wird Integer-Division immer abgeschnitten, so dass 9/2 zu 4 führt. In Python 3 wäre es besser, den Operator // für Ganzzahl-Divisionen zu verwenden, dh (k*k)/3 sollte sein (k*k)//3.

Edit: Sie möchten möglicherweise den Algorithmus, über den Sie gefragt haben, zuordnen: Fastest way to list all primes below N.

+0

Vielen Dank! Das macht viel mehr Sinn. Das bitweise oder hat mich besonders verwirrt. Als 6n + 1 und 6n + 5 auszudrücken, hat es geklickt. Ich hatte gelernt zu beweisen, dass es in der Logik unendlich viele Primzahlen gibt, und auf diese wurde definitiv angespielt. Oh, und dieser Link ist, wo ich diesen Code habe! – Josh