Primzahlenlücken - Optimierungspotential

[nezzcarth]

Hallo
Ich bin seit einer Weile dabei, mich in Python (3.1 und 3.2) einzuarbeiten. Vor ein paar Tagen habe ich versucht mit den Möglichkeiten, die Python so bietet einen möglichst kurzen Primzahlentest (auf Basis einfacher Modulooperationen, nichts elaboriertes wie der Miller-Rabin-Test oder so) zu basteln und bin bei einem Einzeiler gelandet.
Da kam mir die Idee, dass ich darauf aufbauend systematisch Primzahlenlücken (http://oeis.org/A001223) berechnen könnte.

Nach einigem Basteln landete ich irgendwann bei Folgendem:

Code: Alles auswählen

from functools import partial
from operator import mod
from itertools import tee
from time import time

def create_prime_tester(seed=[2]):
    prime_list = seed
    def check_prime(inpt_numb):
        if not inpt_numb in prime_list:
            test_res = all(map(partial(mod, inpt_numb), prime_list))
            if test_res == True:
                prime_list.append(inpt_numb)
            return(test_res)
        elif inpt_numb in prime_list:
            return(True)
    return(check_prime)
        
if __name__ == '__main__':
    is_prime = create_prime_tester()
    inpt = int(input("N: "))
    start_time = time()
    prime_gen = (i for i in range(3, inpt, 2) if is_prime(i) == True)
###
    prime_gen_a, prime_gen_b = tee(prime_gen, 2)
    next(prime_gen_b)
    deltas = ((i, j, i-j) for i, j  in zip(prime_gen_b, prime_gen_a))
    for i, j, k in deltas:
        print ("{0}\t\t{1}\t\t{2}".format(i, j, k)) 
###
    stop_time = time()
    print("Duration:", stop_time - start_time)

(Kommentare und Docstrings habe ich zugunsten der Länge entfernt).
Für etwas größere Zahlen wird das aber schnell relativ langsam (für Primzahlen bis 1.000.000 je nach Rechner 20 1,5 - 2 Std. wenn ich das richtig im Kopf habe).
Teilweise habe ich auch andere Implementationen versucht, konnte aber aufgrund der langen Laufzeit für größere Zahlen oftmals nicht feststellen, ob die Unterschiede natürliche Schwankungen sind, oder nicht.
Für die Berechnung der Differenz hatte ich noch folgenden Kandidaten, der mit einem Generator weniger auskommt:

Code: Alles auswählen

[...]
    first = next(prime_gen)
    while True:
        try:
            second = next(prime_gen)
        except StopIteration:
            break
        delta = second - first
        print ("{0}\t\t{1}\t\t{2}".format(second, first, delta))
        first = second
[...]

(Anstelle des im ersten Source mit ### markierten Abschnittes)

Die konkrete Programmieraufgabe selbst ist nebensächlich. Mir geht es eher um Folgendes:
Wie erreiche ich die optimale Implementierung, unter Beibehaltung der prinzipiellen Idee (Also Primzahlentest per Modulo, und anschließende Differenzbildung und ohne andere Ansätze wie andere Primzahlentest, oder Verwendung mehrerer Threads etc.)?
Kann man das so machen, wie ich das mir gedacht habe, oder sind da Sachen drin, die schlechter Programmierstil, pythonfremd, oder ineffizient sind?
Mir geht es darum, Python anständig zu lernen, und mir keine komischen Dinge anzugewöhnen.

Bin für jede Art von Anmerkung oder Kritik dankbar.

[joshi]

Hallo nezzcarth,

zum Programmierstil kann ich erstmal nur sagen, dass die Klammern nach `return` unnötig sind und du sie weglassen solltest. `return` ist schließlich keine Funktion.

Desweiteren ist deine Methode, die Primzahlen zu erzeugen, ziemlich ineffizient geschrieben. Du speicherst jede Primzahl in einer Liste ab, obwohl du sie später nicht mehr benötigst.
Was du brauchst, ist ja ein generator, der dir nacheinander primzahlen bis zu einem bestimmten Wert liefert. Dazu ist es hilfreich, die Funktion, die Primzahlen überprüft, auszugliedern. Ich hätte da an sowas gedacht (ungetestet!):

Code: Alles auswählen

import math

def is_prime(number):
    if number == 2:
        return True
    for i in range(2, int(math.sqrt(max))+1):
        if number % i == 0:
            return False
    return True

def prime_gen(max):
    yield 2    # 2 ist ein kleiner sonderfall
    for i in range(3, max+1, 2):    # man muss ja nur jede ungerade zahl testen
        if is_prime(i):
            yield i
    raise StopIteration

Für die Berechnung der Lücken ist meines Erachtens nach nur ein Generator nötig, allerdings solltest du die Möglichkeiten der `for`-Schleife nutzen. (wieder ungetestet)

Code: Alles auswählen

first = 2
for prime in prime_gen(max):
    delta = first - prime
    if delta != 0:        # ist nötig, da prime_gen als erste primzahl auch 2 ausspuckt
        #<print delta>
    first = prime

Ich muss zugeben, das mit dem if delta != 0 ist jetzt ein wenig unschön, aber vielleicht fällt noch eine bessere Methode ein

[BlackJack]

@nezzcarth: Also an der Namensgebung kannst Du arbeiten. Zum einen keine Abkürzungen verwenden, also `input_number` statt `inpt_numb`. Das sind ja nur drei gesparte Zeichen, dafür muss man sich aber merken *wie* abgekürzt wurde, jedes mal wenn man den Namen schreibt, und dann kann so etwas auch schnell zu Verwirrungen führen wenn jemand etwas anderes in so einer Abkürzung liest, als gemeint ist. 'numb' ist ja zum Beispiel ein gültiges Wort mit einer Bedeutung. Wird hier jetzt wahrscheinlich nicht verwechselt werden, aber es gibt halt Fälle wo es nicht so eindeutig ist. Das zweite bei Namen: Keine Datentypen in den Namen schreiben. Wenn der sich mal ändern sollte, dann muss man entweder alle Vorkommen von dem Namen anpassen, oder es sind irreführende Namen im Programm. Könnte hier zum Beispiel `prime_list` betreffen, denn da könntest Du mal überlegen ob da nicht ein `set` anstelle einer `list` effizienter wäre. Wenn man das Objekt einfach `primes` nennt, kann man so eine Typänderung vornehmen, ohne dass man den Namen anpassen muss.

``return`` ist keine Funktion -- das sollte man syntaktisch nicht so schreiben als wäre es eine.

Wenn man einen literalen Wahrheitswert aufgrund einer Bedingung zurück gibt, kann man auch gleich das Ergebnis der Bedingung zurück geben. Also bei `check_prime()` kannst Du das ``elif`` weglassen und als letzte Zeile ein ``return input_number in primes`` schreiben.

Explizite Tests auf `True`/`False` sind IMHO schlechter Stil. Bei ``is_prime(i) == True`` wird der erste Teil entweder zu `True` oder zu `False` ausgewertet, also entweder ``True == True`` oder ``False == True``. Und bei den beiden Ausdrücken kommt immer der erste wieder heraus also kann man auf den unnötigen Vergleich mit `True` auch verzichten. Wenn man gegen `False` testen möchte, sollte man statt ``condition == False`` oder ``condition != True`` einfach ``not condition`` schreiben.

Default-Werte werden nur *ein mal* ausgewertet wenn die ``def``-Anweisung ausgeführt wird, und dann wird bei jedem Aufruf das selbe Objekt an den Parameternamen gebunden. Deshalb muss man bei veränderbaren Objekten wie zum Beispiel Listen vorsichtig sein -- wenn das Objekt verändert wird, sind diese Veränderungen bei weiteren Aufrufen noch da! Das heisst die `seed`-Liste enthält alle Primzahlen, die da mal rein gesteckt wurden. Ist in diesem Fall als Cache vielleicht sogar erwünscht. Man sollte sich dessen aber bewusst sein. Und es am besten Dokumentieren, damit andere sehen, dass das Absicht und kein Versehen ist. Das könnte man hier sogar dazu verwenden das Closure weg zu lassen.

Bei der Alternative mit einem Generator weniger könntest Du aus der ``while``-Schleife eine ``for``-Schleife machen, oder!?

Was mir an dem Ansatz nicht so ganz gefällt ist die Abhängigkeit von `check_prime()` davon, dass man es mit den richtigen Zahlen in der richtigen Reihenfolge aufrufen muss. Wenn man da von aussen sowieso zwingend alle ganzen Zahlen (oder zumindest die Primzahlen) in aufsteigender Reihenfolge rein füttern muss, das aber falsch machen könnte, sollte die Generierung der Zahlen auch in der Funktion selber passieren. Man könnte eine Generatorfunktion daraus machen, welche die Primzahlen liefert.

@joshi: Die Primzahlen in der Liste werden sehr wohl benötigt und der Ansatz ist *deutlich* effizienter als Dein Vorschlag.

[numerix]

@nezzcarth: Du müsstest zunächst einmal klären, was genau du erreichen willst, denn davon hängt letztlich ab, ob du mit einem Siebverfahren herangehen solltest oder mit einem probabilistischen Primzahltest. Dass dein Algorithmus schon für n=10⁶ rund 2 Stunden für die Berechnung braucht (wenn ich das richtig verstanden habe), zeigt, dass der eingeschlagene Weg nicht der richtig ist.

Bis zu einer Größenordnung von ca. 10⁸ bist du mit einem vernünftig implementierten Sieb gut bedient. Auf meinem Rechner dauert es deutlich weniger als 10 Sekunden, alle Primzahllücken bis 10⁸ zu ermitteln (mit leichten Geschwindigkeitsvorteilen von Python 3.2 gegenüber Python 3.1). Von 10⁸ bis ca. 10¹⁵ kannst du Teilintervalle in der Größenordnung von 10⁷ bis 10⁸ ebenfalls mit entsprechend angepassten Siebverfahren in wenigen Sekunden durchforsten.

Jenseits von ca. 10¹⁵ wird die Luft dann dünn. Da wäre dann m.E. ein probabilistischer Test geeigneter, allerdings müsste man ihn schon einigermaßen effizient implementieren.

[BlackJack]

@nezzcarth: Eine recht kompakte Implementierung von Deinem Ansatz für einen Iterator der "unendlich" viele Primzahlen generiert, sieht übrigens so aus (Python 2.x):

Code: Alles auswählen

from functools import partial
from itertools import count, ifilter, islice, takewhile


def iter_primes():
    yield 2
    primes = (n * 2 + 1 for n in count(1))
    while True:
        prime = primes.next()
        yield prime
        primes = ifilter(lambda n, p=prime: n % p, primes)


def main():
    result = list(takewhile(lambda n: n <= 1000000, iter_primes()))
    print len(result)


if __name__ == '__main__':
    main()

Der braucht bei mir allerdings auch ca. 20 Minuten.

[derdon]

Kleine unnötige Bemerkung: weder islice noch partial werden in dem Code verwendet (obwohl sie importiert werden).

[nezzcarth]

Hallo,

vielen Dank für die hilfreichen Hinweise und ausführlichen Ratschläge.

@Joshi: Ja, die Sache mit den Klammern bei return sehe ich sofort ein; wollte es wohl zu richtig machen .

Die Primzahlenliste habe ich erst relativ spät eingebaut. Sie stellt meinen Versuch dar, das hier umzusetzen: http://en.wikipedia.org/wiki/Memoization
BlackJacks Anmerkung, dass das effizienter sei, deckt sich mit meiner Erfahrung, dass das Tempo der Berechnung damit ordentlich angezogen hat. Danke trotzdem für deinen Vorschlag.

@BlackJack: In den aller meisten Punkten stimme ich dir sofort zu, und habe entsprechende Überarbeitungen vorgenommen.

Die Sache mit dem eigentlich unnützen elif rührt daher, dass ich mir vor einiger Zeit angewöhnt habe, in Fallunterscheidungen alle Umstände, die mir einfallen explizit abzufragen, alle anderen als undefiniert oder Fehler zu interpretieren und else (das hier zwar nicht syntaktisch ausgedrückt, aber implizit ja irgendwie vorhanden ist) nur zu verwenden, wenn ich mir sehr sicher bin, dass es passt. In diesem Fall hat du aber wahrscheinlich recht, dass ich mir das sparen kann.

Was die Abhängigkeit von "check_prime" bzgl. der Reihenfolge angeht, hast du recht. Das werde ich noch mal in Angriff nehmen. Ich vermute, das dein Code Beispiel (danke auch noch mal dafür) das berücksichtigt. Weil ich allerdings erst mal selbst versuchen möchte, ob ich das hin bekomme, habe ich es mir noch nicht genauer angesehen werde das aber mit Sicherheit noch tun.

Letztendlich bin ich bei Folgendem gelandet:

Code: Alles auswählen

from functools import partial
from operator import mod
from time import time

def check_prime(input_number, prime_cache=[2]):
    if not input_number in prime_cache:
        test_result = all(map(partial(mod, input_number), prime_cache))
        if test_result:
            prime_cache.append(input_number)
        return test_result
    return input_number in prime_cache
        
if __name__ == '__main__':
    limit = int(input("N: "))
    start_time = time()
    prime_generator = (i for i in range(3, limit, 2) if check_prime(i))
    first = next(prime_generator)
    for second in prime_generator:
        delta = second - first
        print ("{0}\t\t{1}\t\t{2}".format(second, first, delta))
        first = second
    stop_time = time()
    print("Duration:", stop_time - start_time)

Schon wesentlich kürzer. Seltsamer Weise hatte ich das Gefühl, dass das Einsparen der Closure etwas Geschwindigkeit gebracht hat, die Umstellung auf einen Generator das aber wieder ausgeglichen hat. Werde dem noch mal nachgehen

@Numerix: Mir ging es in erster Linie darum, mal zu sehen, wie ich eine konkrete Idee in Python sauber und effizient umsetzen kann. Einen wirklichen Nutzen in dem Sinne, dass ich das Programm für zum Beispiel wissenschaftliche Zwecke oder so bräuchte, besteht nicht, sondern nur in dem Übungseffekt für mich (mal abgesehen davon, dass ich Primzahlen wie viele andere spannend finde und einfach neugierig bin). Was die Laufzeit anbelangt ist mir erst nachträglich aufgefallen, dass ich euch einige wesentliche Faktoren unterschlagen habe.
Allerdings beträgt auch die Laufzeit unter weniger besonderen Bedingungen mit c.a. 55 Minuten noch wesentlich mehr, als bei dir. Es stimmt also wohl, dass da grundsätzlich etwas falsch läuft. Wenn ich dich richtig verstehe, siehst du das wesentliche Problem im Primzahlentest. Ich werde in der nächsten Zeit mal schauen, dass ich einen effizienteren Test implementiere. Das setzt nur voraus, dass ich mich ggf. vorher in die dahinter stehende Theorie einarbeiten muss.

Ein Lerneffekt hat sich dennoch schon eingestellt, und eure Tipps haben mir in jedem Fall was gebracht

EDIT: Kleine Codeänderung vorgenommen.

[BlackJack]

@nezzcarth: Ich denke numerix wollte, zumindest für den Bereich bis 1 Mio., nicht unbedingt auf einen besseren Primzahltest für einzelne Zahlen hinaus, sondern dass man gar nicht testet sondern die Zahlen mit einem Siebverfahren bestimmt.

Ansonsten hat das IMHO nur begrenzten praktischen Wert in Python. Das fällt IMHO unter "number crunching" und ist nicht der bevorzugte Einsatzbereich für eine dynamische, interpretierte Sprache. Zumindest nicht wenn man es komplett in der Sprache umsetzt. Wenn man bei so etwas Geschwindigkeit heraus holen will, sollte man auf externe Bibliotheken zurück greifen. Siebverfahren sollte man mit `numpy`-Arrays zu Beispiel beschleunigen können. Oder es in C implementieren und per `ctypes` aufrufen.

[numerix]

Wenn ich dich richtig verstehe, siehst du das wesentliche Problem im Primzahlentest.

Wenn ich jetzt "ja" sage, könntest du das missverstehen. Ob dein Testverfahren "gut" oder "schlecht" ist, hängt zunächst einmal von deinem Ziel ab. Geht es darum, einzelne (vorzugsweise ziemlich große) Zahlen auf Primalität zu prüfen, dann würde man eher mit einem probilistischen Test wie dem von dir schon genannten Miller-Rabin-Test an die Sache herangehen. Den zu implementieren ist nicht sooo schwierig und wenn du Spaß an Primzahlen hast, kann es durchaus reizvoll sein, eine möglichst performante Implementation zu erreichen.

Wie gut dir das gelungen ist, könntest du dann z.B. hier überprüfen:
https://www.spoj.pl/problems/PRIC/
https://www.spoj.pl/problems/PON/
https://www.spoj.pl/problems/PAGAIN/ (mit Python nur schwer schaffbar)

Wenn es hingegen darum geht, eine größere Menge an aufeinanderfolgenden Primzahlen zu berechnen - und genau das brauchst du ja, um (systematisch) Primzahllücken zu ermitteln - dann wäre ein Siebverfahren der richtige Weg. Genial ist das Sieb des Eratosthenes: Einfach zu verstehen, einfach zu implementieren und dennoch enorm effizient (wenn man es vernünftig macht). Auch dazu einige Aufgaben (die letzten beiden sind nur mit psyco und reichlich Optimierung schaffbar.):

https://www.spoj.pl/problems/PRIME1/
https://www.spoj.pl/problems/TPRIMPER/
https://www.spoj.pl/problems/CPRIME/
https://www.spoj.pl/problems/TDPRIMES/