Nach einer etwas längeren Pause geht es heute mit der Projecteuler-Serie weiter. In diesem Artikel stelle ich meine Lösung des dritten Problems dar. Wie bereits bei den vorangehenden Teilen gibt es auch hier wieder alternative Ansätze.

Problemvorstellung

Das dritte Problem von Projecteuler.net dreht sich allgemein gesprochen um die Primfaktorzerlegung. Aus grauer Schulzeit wissen wir, dass jede natürliche Zahl in Primfaktoren zerlegt werden kann. Des Weiteren haben Primzahlen die berüchtigte Eigenschaft nur durch sich selber und durch die Eins teilbar zu sein. Danach hat jede natürliche Zahl mindestens zwei Faktoren, wobei die Eins als trivialer Faktor beinahe vernachlässigt werden kann. Nicht-Primzahlen, also die zusammengesetzten Zahlen, können folglich in ihre kleinsten Bestandteile, die Primzahlen getrennt werden.

Ich muss zugeben, dass es nun wirklich keine wissenschaftliche Erklärung und kein Beweis für die Existenz einer exakten Primfaktorzerlegung ist, was ich aufgeführt habe. Den Zweck einer netten Einführung erfüllt es jedoch und ich bin glücklich . Sollte jemand weiterführende Informationen benötigen, so möge Er oder Sie sich im Web umschauen oder einen Hinweis in den Kommentaren hinterlassen, dann lässt sich sicher noch etwas machen.

Also wir wissen nun, wohin das dritte Problem führen wird, doch wie lautet die Aufgabenstellung genauer? Die Auflösung:

Finde den größten Primfaktor der Zahl 600851475143.

Ja, das war es schon! Hört sich doch ganz machbar an, oder? Das haben sich 56838 der derzeit 100370 Hobbymathematiker auch gedacht und die Aufgabe erfolgreich gelöst. Im Vergleich zum ersten und zweiten Problem entspricht dies schon deutlich weniger eingereichten Lösungen. Bisher sind es grob abgeschätzt etwa 50% aller Projecteuler.net Aktivisten. Schlussendlich kann diese Entwicklung auf einen steigenden Schwierigkeitsgrad oder einen sinkenden Motivationsgrad hindeuten, also ein typisches Glas-halb-leer bzw. Glas-halb-voll Problem, welches auf in einem anderen Artikel gelöst wird .

Problemeingrenzung

Kommen wir nun zum Kern der Geschichte: Das erste Durchlesen dieses Problems hinterließ bei mir zwei Gedanken. Der erste Gedanke war, wie finde ich raus, ob eine Zahl nun eine Primzahl ist oder nicht und der zweite Gedanke beschäftigte sich mit dem Prüfen des größten Teilers.

Das erste Teilproblem gehe ich im folgenden mit dem Sieb des Erathosthenes (Wikipedia) an.

Das zweite Teilproblem wird ganz simpel gelöst, indem die gegebene Zahl auf Teilbarkeit durch die gefundenen Primzahlen untersucht wird.

Problemlösung

Das Sieb des Erathostenes ist eine recht einfache Methode alle Primzahlen bis zu einer bestimmten Zahl N zu finden. Wie im verlinkten Wiki-Artikel schön dargestellt ist, schreibt man alle Zahlen bis einschließlich N auf und beginnt das Sieben. Dazu fängt man mit der kleinsten bekannten Primzahl, also der zwei, an und streicht alle Vielfachen dieser Zahl, dann macht man dasselbe mit der nächsten nicht gestrichenen Zahl nach der Zwei weiter. Auf diese Weise bleiben am Ende alle Primzahlen bis N übrig.

Der langen Rede kurzer Sinn ist in Python-Code im Folgenden dargestellt.

import time as t
import numpy as np
from math import ceil, sqrt

def eratosthenes_direct(N):
    """
    Berechnet alle Primzahlen bis N mit dem Sieb des Eratosthenes
    """
    primes = np.arange(3,N+1,2)  # only uneven figures
    primes = np.insert(primes,0,2,0) # ... and two

    testFigs = np.arange(3,np.ceil(np.sqrt(N))+2,2)

    for i in testFigs:
        if (primes[((i+1)/2-1)] != 0):
            noPrime = i*i
            while noPrime <= N:
                primes[((noPrime+1)/2-1)] = 0
                noPrime = noPrime + 2 * i

    return np.trim_zeros(np.sort(primes))

def eratosthenes_vect(N):
    primes = np.arange(3,N+1,2)  # only uneven figures
    primes = np.insert(primes,0,2,0) # ... and two

    nUneven = np.size(primes) - 1.0

    for i in np.arange(3,np.ceil(np.sqrt(N))+2,2):
        primes[(i*i+1)/2-1:nUneven+1.0:i] = 0;

    return np.trim_zeros(np.sort(primes))

In dem Code-Abschnitt habe ich gleich zwei alternative Beispiele für den Sieb des Eratosthenes aufgeführt. Beide Funktionen gehen nach dem gleichen Prinzip vor; der kleine Unterschied ist beim Streichen der Vielfachen zu finden.

Da von vornherein klar ist, dass keine gerade Zahl außer der Zwei Primzahl sein kann, werden nur die ungeraden Zahlen bis N im Vektor primes gespeichert. Bevor ich den Code getippt habe, führte ich das Sieben auf dem Papier aus und bin dann darauf gekommen, dass die größte, auf Vielfache zu prüfende Primzahl die Wurzel von N sein muss. Alle Zahlen die größer als die Wurzel von N sind, müssen demnach Primzahlen sein oder werden beim Sieben gestrichen.

Die an zweiter Stelle aufgeführte Funktion nutzt die Möglichkeiten von Numpy aus, indem keine while-Schleife verwendet wird, sondern das Streichen über eine Indexauswahl des Vektors stattfindet.

Zuletzt werden mit den Numpy-Funktionen sort und trim_zeros die Primzahlen von den Nullen getrennt.

So weit so gut! Nun liefert uns das Sieb einen Vektor an Primzahlen bis zu einer bestimmten Zahl N. Um die Projecteuler-Aufgabe zu lösen, muss noch der größte Primfaktor gefunden werden:

def solver(N,primes):
    rebuild = 1.0
    prob    = N
    for i in primes:
        while (N%i == 0):
            N = N/i
            found = i
            rebuild = rebuild * found

            if (rebuild/prob > 0.9):
                print(rebuild)
                return found

    return 0.0

Der Solver durchläuft alle Primzahlen, prüft Sie auf Teilbarkeit und, wenn die Zahl ein Teiler ist, so wird sie als ein Lösungskandidat gespeichert. Die while-Schleife soll doppelte Primzahlen ebenfalls erfassen.

Der Solver wurde von mir nicht weiter validiert und könnte durchaus noch das eine oder andere Problem bereiten, aber für die von mir getesteten Zahlen hat alles funktioniert.

Lösungsanalyse

Kommen wir nun zur Anwendung des Siebs und des Solvers.

Zunächst muss eine Grenze für das Sieb festgelegt werden. Da ich den Solver so geschrieben habe, dass er mir nur einen Lösungskandidaten ausgibt, wenn die Zahl von Projecteuler rekonstruiert werden konnte, habe ich als Grenze zunächst 1000 ausprobiert. Das war zu wenig und ich erhöhte die Grenze auf 10000, was eine Lösung ausspuckte.

Der Aufruf:

N  = 10000
Pe = 600851475143
#
t0 = t.clock()
primes = eratosthenes_direct(N)
prim = solver(Pe,primes)
tn = t.clock()

print("Zeit fuer Primzahlensieb: %.15f" %(tn-t0))
print("Primzahl: ")
print(prim)

t0 = t.clock()
primes2 = eratosthenes_vect(N)
prim2 = solver(Pe,primes2)
tn = t.clock()

print("Zeit fuer Primzahlensieb: %.15f" %(tn-t0))
print("Primzahl: ")
print(prim2)

Und die Ausgabe:

600851475143.0
Zeit fuer Primzahlensieb: 0.080000000000000
Primzahl:
6857
600851475143.0
Zeit fuer Primzahlensieb: 0.030000000000000
Primzahlen:
6857

Abgesehen davon, dass das Ergebnis die richtige Lösung für das dritte Projecteuler-Problem ist, kann festgestellt werden, dass der Sieb ohne die while-Schleife um mehr als das doppelte schneller ist. Das war zu erwarten, denn die Numpy-Funktionen werden sicherlich nicht von Anfängern, wie ich einer bin, geschrieben.

Fazit

Damit wäre das dritte Problem gelöst. Auch diesmal gibt es eine kleine Neuerung im Blog: mit dem Sieb des Eratosthenes wurde der erste Algorithmus vorgestellt. Ich denke, dass ich in nächster Zeit wieder etwas mehr Zeit für den Blog haben werde und daher öfter etwas schreiben kann.

Jetzt sind aber erstmal die Leser wieder an der Reihe mit Verbesserungsvorschlägen und weiteren Lösungsansätzen. Ich freue mich darauf.

Artikel aus der selben Kategorie:

• » Projecteuler: Lösung zu Problem 10
• » Projecteuler: Lösungen in Serie
• » Projecteuler: Lösung zu Problem 4
• » Projecteuler: Lösung zu Problem 7
• » Projecteuler: Lösung zu Problem 5