Seite 3: Die Suche nach den Punkten

Eine Suche nach solchen Punkten (x,y) mit x bis zu 100.000000 (Hundertmillionen) dauert auf der gleichen Hardware wie bei den vorherigen Beispielen auch nur 5 Sekunden. Jedoch liefert die Routine neben den zwei richtigen Lösungen auch ein paar, die gar keine sind:

• richtig:x=0 y=2
• richtig: x=1, y=1
• falsch: x=65536, y=4294967294
• falsch: x=83913453, y=824417755

Die NumPy-Datentypen und das Numba-Paket sind aufeinander abgestimmt. Programmierer, die von der Optimierung profitieren möchten, sollten deshalb NumPy-Datentypen verwenden. Der größtmögliche NumPy-Integer ist uint64 (unsignierter 64 Bit Integer). Er kann Werte annehmen, die bis zu 2⁶⁴-1 =18.446.744.073.709.551.615 groß sind. Da jedoch die Variable x eine Schleife von 0 bis 100.000000 durchläuft und es gemäß der Kurvengleichung hoch 6 gerechnet wird, sprengt es die Grenze. Es entsteht ein Überlauf, der falsche Zahlen und falsche Ergebnisse produziert. Leider ist es für Entwickler alles andere als einfach, damit umzugehen. Unter Verwendung des Pakets gmpy2 für hohe Zahlenpräzision tauchen solche falschen Lösungen zwar nicht mehr auf, jedoch dauert die in Listing 6 implementierte Routine statt 5 Sekunden nun 842 Sekunden (14 Minuten).

Nimmt man die JIT-Optimierung aus Listung 6 heraus, dauert sie etwa gleich lang, nämlich 769 Sekunden (13 Minuten). Am schlechtesten ist die Performance, wenn die Optimierung aus Listing 5 entfernt wird, man also nichtoptimiert die NumPy-Datentypen anstelle der gmpy2-Datentypen verwendet. In diesem Fall treten keine Falschlösungen auf, allerdings dauert die Routine statt 5 Sekunden nun 6667 Sekunden (1,8 Stunden). Es zeigt sich, dass die hauptsächliche Optimierung in gmpy2 erfolgt. Mithilfe des in Python verfügbaren multiprocessing-Pakets ist es möglich, die Routine wieder stark zu beschleunigen – nämlich runter auf 49 Sekunden (siehe Funktion find_integer_solutions_gmpy2_multiprocessing im Github-Repository ).

@jit('void(uint64)')
def findIntegerSolutionsGmpy2(limit: np.uint64):
    for x in np.arange(0, limit+1, dtype=np.uint64):
        x_ = mpz(int(x))**2
        y = x_**3-mpz(4)*x_+mpz(4)
        if gmpy2.is_square(y):
            print([x,gmpy2.sqrt(y),y])

Listing 6: Die Routine aus Listing 5 auf gmpy2-Päzisionsganzzahlen umgestellt.

Alternativ kann auch das Werkzeug Mathematica (Wolfram) Abhilfe schaffen. Es ist dazu in der Lage, solche Berechnungen mit hoher Präzision durchzuführen (Listing 7). Allerdings dauert die Berechnung in diesem Fall 13856 Sekunden (3,8 Stunden).

F[x_] := x^6 - 4 x^2 + 4;
sQ[n_] := FractionalPart@Sqrt[n + 0``1] == 0;
For[i = 0, i < 1000000000, 
 i++, {x = i; y = F[x]; 
  If[sQ[y] == True, Print[{x, Sqrt[y]}], Continue]}]

Listing 7: Die Routine aus Listing 5 in Mathematica liefert nur korrekte Lösungen.

Schließlich soll es auch noch um die Implementierung von CUDA gehen: In diesem Fall läuft der Code auf einer GPU, wozu Entwickler den Code aus Listing 5 entsprechend anpassen müssen. Listing 8 zeigt das Ergebnis der Anpassung. Hierbei kommt für die optimierte Funktion statt des @jit-Dekorators der @cuda.jit-Dekorator zum Einsatz. Wenn das Programm die Funktion aufruft, lassen sich die Blockgröße 128 und die Anzahl der Threads pro Block 255 übergeben. Die beiden Werte ergeben miteinander multipliziert die Gesamtzahl der Threads. Die Blockgröße legt fest, wie viele Threads sich einen bestimmten Bereich des gemeinsamen Speichers teilen. Den Code in Listing 8 hat der Autor des Artikels auf einer Maschine mit 128 GByte RAM, einem Intel Xeon CPU E5-2630 v4, 2.20GHz und zwei Grafikkarten vom Typ Tesla V100 mit jeweils 16 GByte RAM ausgeführt. Die Ausführungsdauer liegt bei 316 Sekunden (5 Minuten). Die Performance ist beeindruckend, allerdings enthält das Ergebnis erneut zwei falsche Lösungen.

%%time
from numba import jit, cuda
import numpy as np
from math import sqrt

@cuda.jit
def findIntegerSolutionsCuda(arr):
    i=0
    for x in range(0, 1000000000+1):
        y = float(x**6-4*x**2+4)
        sqr = int(sqrt(y))
        if sqr*sqr == int(y):
            arr[i][0]=x
            arr[i][1]=sqr
            arr[i][2]=y
            i+=1

arr=np.zeros((10,3))
findIntegerSolutionsCuda[128, 255](arr)

print(arr)

Listing 8: Für CUDA umgeschriebene Routine

Um die Auswirkung der Optimierungen hinsichtlich Performance versus Präzision transparent zu machen, stellt Tabelle 2 die Werte gegenüber. Die Auswirkung von Präzision auf Performance ist immens (Faktor 168 bis 1333).

Tabelle 2: Vergleich zu den Auswirkungen der Optimierung.
*Ausführung auf Surface Pro 6 (16GB RAM, Intel i7-8650U, 1.90GHz, 4 Kerne, 8 logische, Intel UHD Graphics 620)