Python: Richtige JIT-Optimierung – damit nichts schiefgeht

07.06.2022 10:32 Uhr Eldar Sultanow

JIT-optimierter Code beschleunigt Python-Programme auf GPUs rasant. Doch wer nicht aufpasst, zerstört schnell den Performancegewinn oder rechnet sogar falsch.

Forscherinnen und Entwickler sehen sich heute oftmals mit Fragestellungen konfrontiert, die immer komplexer werden – auch die Datenmengen nehmen an Komplexität zu. Die Leistung der Hardware wird aber ebenfalls immer besser und steht sogar per Knopfdruck in der Cloud zur Verfügung [1]. Allerdings ist es gar nicht so einfach, diese Leistung optimal auszuschöpfen. Nicht selten laufen Berechnungsroutinen auf leistungsstarker Hardware, nutzen dabei aber nur einen Bruchteil dieser Leistung. Wer on Demand beispielsweise in der AWS-Cloud die EC2 P3-Instanz [2]vom Typ p3.8xlarge (diese beinhaltet vier GPUs) verwendet, bezahlt pro Stunde 12,24 US-Dollar.

Wie sträflich wäre es, dort eine Berechnungsroutine laufen zu lassen, die weder die JIT-Optimierung noch die GPU nutzt: Das dauert unnötig länger und kostet dementsprechend mehr Geld. Von welchen Faktoren hier die Rede ist, wie sehr viel schneller eine JIT- oder GPU-optimierte Routinen läuft und wie Entwicklerinnen und Entwickler sie mit Python umsetzen können, zeigt dieser Artikel. Ausgewählte Beispiele veranschaulichen die Fallstricke, die auftreten können.

Das nötige Handwerkszeug

Die Programmiersprachen C++ und Python sind sehr gut geeignet, um rechenintensive, leistungsoptimierte Routinen zu schreiben – erst im letzten Jahr hat Nvidia Neuerungen für beide Programmiersprachen eingeführt [3]. Dieser Artikel konzentriert sich auf Python, betrachtet dabei aber auch Alternativen über den Tellerrand hinaus (Mathematica, C++). Daher benötigen Entwicklerinnen das entsprechende Handwerkszeug:

Eine Entwicklungsumgebung (IDE): Die hier vorgestellten Beispiele, Programme und Codefragmente wurden in der IDE Visual Studio Code entwickelt. Es gibt auch weitere gut geeignete Umgebungen wie PyCharm.
Eine Python-Installation: Die hier verwendete Installation umfasst die im November letzten Jahres veröffentlichte Anaconda-Version Anaconda Individual Edition 2021.11, die Python 3.9.7 mitbringt.
Entsprechende Pakete: zum Ausprobieren der hier entwickelten Routinen müssen die Pakete Numba (JIT-Compiler zur Übersetzung von Python-Code in schnellen Maschinencode), CUDA Toolkit (GPU-beschleunigte Bibliotheken) Pyculib (Python bindings für CUDA), math (Paket mit mathematischen Funktionen), pandas (Paket zur Verarbeitung und Speicherung/Export von Daten), NumPy (Paket für numerische Berechnungen), multiprocessing (Paket für Prozess-basierte Parallelisierung) und gmpy2 (Paket für hoch performante Multipräzisionsarithmetik) vorhanden sein; time und sys für Zeitmessung und Systemfunktionen sind standardmäßig in Python vorhanden.
Codeassistent und Qualitätsprüfer: Pylint erkennt Probleme, Fehler im Code sowie Verstöße gegen die Konventionen und Empfehlungen der Python-Community und weist auf diese hin.

Anschauungsmaterial: Mengolis Sechs-Quadrate-Problem

Angeregt durch den französischen Mathematiker Jacques Ozanam beschäftigte sich sein italienischer Kollege Pietro Mengoli in den 1670er Jahren mit speziellen diophantischen Gleichungen – sie lassen lediglich ganzzahlige Lösungen zu. Ein bekanntes Problem, das von solchen Gleichungen handelt, ist das Sechs-Quadrate-Problem [4]: Finde drei natürliche Zahlen x, y, z, deren Differenzen z-y, z-x, y-x jeweils Quadratzahlen sind, wobei die Differenzen z²-y², z²-x², y²-x² der Quadrate ebenfalls Quadratzahlen sind.

Die Ergebnisse umfassen sehr große Zahlen. Ozanam hatte damals schon eine Lösung präsentiert, nämlich x=1873432, y=2288168 und z=2399057. Derartige Probleme in höheren Zahlenregionen sind als Anschauungsmaterial ideal, da es Fallstricke in der JIT-optimierten Programmierung mit sehr großen Zahlen gibt, wie sich später noch zeigt.

Die Bedeutung heutiger Rechenleistung

Bemerkenswert ist, dass Ozanam bereits im 17. Jahrhundert derartig großen Zahlen gefunden hat, für die diese Gleichung erfüllt ist. Hatte er doch keine technischen Hilfsmittel wie Computer, Programmiersprachen oder Ähnliches zur Hand. Unsere heutigen Werkzeuge sind in der Lage, solch aufwendige Probleme zu lösen. So hat beispielsweise ein Team am Massachussetts Institute of Technology (MIT) im Jahr 2019 mit x=-80538738812075974, y=80435758145817515 und z=12602123297335631 eine Lösung für die diophantische Gleichung x³+y³+z³=42 gefunden. Unter den Zahlen bis 100 war die 42 die letzte noch offene Zahl, die als Summe dreier Kubikzahlen darstellbar ist. Es wundert nicht, dass das Interesse an diesem Rätsel auch bei Fans des Science-Fiction-Klassikers "Per Anhalter durch die Galaxis" groß war. Das Team hat für die Berechnung die ungenutzte Leistung von mehr als einer halben Million Heim-PCs [5] verwendet. Derartige Gleichungen gehören zu dem seit 160 Jahren ungelösten "Summe dreier Kubikzahlen"-Problem. Es gibt weitere Fragestellungen, die ohne hohe Rechenpower auch in Hunderten von Jahren nicht zu bewältigen wären. Hierzu zählen vor allem Aufgaben aus dem Umfeld der künstlichen Intelligenz (KI) wie Spracherkennung, Bilderkennung oder Mustererkennung in massenweise numerischen Daten.

Algorithmus zur Suche von Quadratzahlen

Zur Vereinfachung liegt der Fokus auf den drei Quadratzahlen x², y² und z², deren Differenz wieder Quadratzahlen sind. Hier gibt es auch Ergebnisse, die nicht ganz so groß sind: x=153, y=185 und z=697 oder x=520, y=533 und z=925. Ein kleiner Quercheck 925²-533² ergibt 571536, was tatsächlich eine Quadratzahl ist, nämlich 756². Nachzuprüfen, ob die anderen Differenzen ebenfalls stimmen, sei an dieser Stelle der Leserschaft überlassen. Eine große Anzahl solcher Treffer bis in die Zahlenregion von 12 Millionen wurden mithilfe einer Just-in-Time-optimierten Routine (JIT) ermittelt und ist in einem GitHub-Repository [6] (siehe Datei pythagorean_12000000.csv) verfügbar. Der dafür zuständige Algorithmus ist auf StackExchange Mathematics [7] beschrieben, für JIT optimiert und ebenfalls im GitHub-Repository (siehe Datei pythagorean_gendata_jit.py) verfügbar.

@jit('void(uint64, uint64[:])')
def generateData(limit: np.uint64, triples: np.ndarray) -> np.uint32:
    A=np.zeros(limit, dtype=np.uint64)
    B=np.zeros(limit, dtype=np.uint64)
    rows = np.uint32(0)
    for w in np.arange(1, limit+1, dtype=np.uint64):
        count = np.uint32(0)
        for a in np.arange(1, w+1, dtype=np.uint64):
            sqr = np.uint64(sqrt(w*w-a*a))
            if sqr*sqr == w*w-a*a:
                count+=1
                A[count]=a
                B[count]=np.uint64(sqrt(w*w-a*a))
        if count>1:
            for i in np.arange(1, count+1, dtype=np.uint64):
                for j in np.arange(1, count+1, dtype=np.uint64):
                    if i!=j:
                        x=np.uint64(A[i])
                        t=np.uint64(B[i])
                        s=np.uint64(A[j])
                        z=np.uint64(B[j])
                        if z > x:
                            y = np.uint64(sqrt(z*z-x*x))
                            if y*y == z*z-x*x:
                                arr = np.array([x, y, z, s, t, w], dtype=np.uint64)
                                triples[rows] = arr
                                rows+=1
    return rows

Listing 1: Optimierter Algorithmus zur Lösung des Sechs-Quadrate-Problems

Wer die Routine in Listing 1 zur Suche von Lösungen startet, die Werte bis zu 50000 annehmen können (limit=50000), benötigt dazu auf einem HP Z-Book (32 GByte RAM, Intel Core i7-9750H, 2.60GHz, 6 Kerne, 12 logische, NVIDIA Quadro T1000) 4,8 Sekunden und findet 1074 Tripel für x, y und z. Führt man die Routine ohne Optierung auf demselben Rechner aus, dann benötigt sie 1082,8 Sekunden (18 Minuten), um diese 1074 Tripel zu ermitteln. Demnach hat die JIT-Optimierung eine Beschleunigung um etwa den Faktor 225 erbracht.

Wie die Optimierung zustande kommt

Wie kommt nun die Optimierung zustande? Dafür ist das Paket Numba zuständig – über from numba import jit lässt sich der dort definierte @jit-Dekorator importieren. Die Schwierigkeit in der JIT-optimierten Programmierung besteht darin, dass die Verwendung von Befehlen eingeschränkt ist und Entwickler nicht einfach beliebigen Code optimieren können. Deshalb muss er optimierbare Operationen in eine dafür gesondert vorgesehene Routine auslagern. Im vorliegenden Fall ist es die Funktion generateData, die mit dem @jit-Dekorator versehen ist.

Die Routine sucht mithilfe dreier verschachtelter for-Schleifen, in denen sie einfache Operationen wie Addition, Subtraktion, Multiplikation und Wurzelziehen verwendet, nach den Zahlen x, y, z. Am Ende gibt sie die Anzahl der gefunden Ergebnisse zurück. Das Speichern der gefundenen Lösungen in eine CSV-Datei erfolgt außerhalb dieser Routine. Fügen Programmierer derartigen, beispielsweise auf pandas basierenden Code zur Erzeugung von CSV-Dateien in die per @jit-Dekorator versehene Methode ein, erhalten sie beim Aufruf des Python-Programms einen Fehler. Der Grund dafür ist, dass derartiger Code nicht optimierbar ist. Allerdings existiert eine Reihe von Fallstricken, die zwar keinen Fehler, aber einen enormen Performance-Verlust und falsche Ergebnisse verursachen.

Fallstrick 1: Array-Dynamiken sind Performance-Killer

Der Code in Listing 1 ist so beschaffen, dass er ein leeres, bereits mit vorgegebener Größe initialisiertes Array an die optimierte Routine übergibt. Das Programm befüllt dieses Array innerhalb der Route selbst. Alternativ könnten Entwickler dem Gedanken verfallen, das Ergebnis-Array dynamisch mit den gefundenen Lösungen wachsen zu lassen. Eine Suche von Lösungen, deren Werte x, y, z maximal 50000 annehmen, ergibt viel weniger Ergebnisse als die festgelegte Grenze von 50000 – in der Tat sind 1074 Treffer deutlich weniger als 50000. Auf diese Weise würde das Programm auch keinen Platz verschwenden.

Der Programm-Code im Listing 2 realisiert diese Platzersparnis und hält sich dabei an die Vorgabe, dass er an die optimierte Routine ein initialisiertes Array übergeben soll. Allerdings zerstört diese Vorgehensweise den gesamten Performancegewinn und eine Suche (wieder mit limit = 50000) dauert nun 1203,3 Sekunden (20 Minuten). Es ist also von Vorteil, wenn sich Programmierer auf diese "Platzverschwendung" einlassen. Wer dieses Beispiel selbst nachstellen möchte, findet den vollständigen Code im GitHub-Repository (siehe Datei pythagorean_gendata_jit_antipattern.py).

arr = np.array([x, y, z, s, t, w], dtype=np.uint64)
old_size = triples.shape
rows = np.uint32(old_size[0])
cols = np.uint32(old_size[1])
triples.resize((rows + 1, cols), refcheck=False)
triples[rows] = arr

Listing 2: Fallstrick 1 – "Teure" Operationen in der optimierten Routine.

Bei der Gegenüberstellung der Laufzeiten, wie sie in Tabelle 1 zu sehen sind, ergibt sich ein beeindruckender Vergleich: Die korrekt umgesetzte Optimierung bringt eine Beschleunigung um den Faktor 225.

Keine Optimierung	Optimierung	Optimierung mit Fallstrick 1
1082,8 Sekunden	4,8 Sekunden	1203,3 Sekunden

Tabelle 1: Gegenüberstellung der Laufzeiten

Fallstrick 2: Dynamische Schleifen verhindern die Optimierung

Setzen Programmierer eine lang laufende Routine ein, deren Laufzeit nicht nur Tage, sondern Wochen dauert, ist es wünschenswert, in gewissen Zeitabständen Zwischenergebnisse zu sichern – beispielsweise in einer CSV-Datei. Der Schreibvorgang darf aber nicht innerhalb der optimierten Routine, sondern muss außerhalb erfolgen. Listing 3 zeigt eine Routine mit drei ineinander verschachtelten for-Schleifen zur alternativen Generierung potenzieller Tupel für das Sechs-Quadrate-Problem. Der vollständige Code ist im GitHub-Repository (siehe Datei pythagorean_gendata2_nofile_jit.py) verfügbar. Dieser Code lässt sich per JIT gut optimieren und bietet die die erwartete Beschleunigung um dreistellige Faktoren.

@jit('void(uint64)')
def generateData(limit: np.uint64):
    for t in np.arange(0, limit+1, dtype=np.uint64):
        for s in np.arange(0, t, dtype=np.uint64):
            for u in np.arange(0, limit+1, dtype=np.uint64)

Listing 3: Ineinander geschachtelte For-Schleifen mit festen Bereichen kann JIT optimieren.

Speichern Entwickler die gefunden Ergebnisse in einer Datei zwischen, entfällt dabei die gesamte Optimierung. Jede Fundstelle wird an den (außerhalb der JIT-Optimierung laufenden) Aufrufer zurückgegeben, der sie wiederum in eine Datei schreibt (Listing 4). Der Mechanismus hebt die gesamte JIT-Optimierung auf. Im GitHub-Repository [8] liegt der vollständige Code zum Ausprobieren bereit (siehe Datei pythagorean_gendata2_jit_antipattern.py).

@jit('void(uint64, uint64[:])')
def generateData(limit: np.uint64, triplet: np.ndarray) -> np.ndarray:
    s0 = triplet[0]
    t0 = triplet[1]
    u0 = triplet[2]
    for t in np.arange(t0, limit+1, dtype=np.uint64):
        for s in np.arange(s0, t, dtype=np.uint64):
            for u in np.arange(u0, limit+1, dtype=np.uint64):
               �
               if sqr*sqr == t_s:
                   return np.array([s, t, u, ss, tt, uu, t_u, t_u_s, t_s], dtype=np.uint64)

Listing 4: Aufheben der Optimierung durch dynamische Schleifen (Fallstrick 2).

Fallstrick 3: Datentypen in der optimierten Routine

Der dritte Fallstrick zeigt sich anhand eines in der Kryptographie verbreiteten Problems: Elliptische Kurven. Eine häufige Fragestellung, die direkt damit im Zusammenhang steht, lautet: Liegen auf der Kurve rationale oder sogar ganzzahlige Punkte? Die Frage wurde beispielsweise auf StackExchange Mathematics [9]für die Kurve y²=x⁶−4x²+4 gestellt und beantwortet. Es sei an dieser Stelle erwähnt, dass es mächtigere mathematische (beispielsweise probabilistische) Methoden für die Suche von rationalen Punkten auf elliptischen Kurven gibt. Zu vereinfachten Illustrationszwecken des Performancegewinns via JIT und GPU bleibt dieses Beispiel beim einfachen Bruteforce.Die Routine in Listing 5, die in diesem Code korrekt optimiert ist, dient zur Suche von ganzzahligen Lösungen dieser Kurve.

@jit('void(uint64)')
def findIntegerSolutions(limit: np.uint64):
    for x in np.arange(0, limit+1, dtype=np.uint64):
        y = np.uint64(x**6-4*x**2+4)
        sqr = np.uint64(np.sqrt(y))
        if np.uint64(sqr*sqr) == y:
            print([x,sqr,y])

Listing 5: Optimierte Routine zur Suche von ganzzahligen Punkten auf der Kurve y²=x⁶−4x²+4

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 [10]).

@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).

Kein JIT, numpy 64 Bit Integer	JIT-Optimierung, numpy 64 Bit Integer	JIT-Optimierung, gmpy2 Präzisions-Ganzzahlen	Kein JIT, gmpy2 Präzisions-Ganzzahlen	Mathematica*	CUDA.JIT	gmpy2, multicore *	C++ mit GMP *
6667 Sekunden	5 Sekunden	842 Sekunden	769 Sekunden	13856 Sekunden	316 Sekunden	49 Sekunden	7 Sekunden
Nur korrekte Lösungen	Korrekte und falsche Lösungen	Nur korrekte Lösungen	Nur korrekte Lösungen	Nur korrekte Lösungen	Korrekte und falsche Lösungen	Nur korrekte Lösungen	Nur korrekte Lösungen

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)

Schnell und präziser: C++

Und welche Ergebnisse bekommen Programmierer, die lieber C++ als Python einsetzen? Dieser Ansatz ist zwar eher unhandlich und sicher nicht so populär wie die Verwendung von Python. Das Ergebnis ist jedoch erstaunlich. Listing 9 zeigt den Code (der komplette Code steht ebenfalls auf GitHub [11]zum Download bereit).

Die Performance der Bruteforce-Suche von ganzzahligen Punkten auf der Kurve y²=x⁶−4x²+4 als Balkendiagramm.

Um Präzisionsverlust zu vermeiden, kam bei diesem Beispiel die Bibliothek GMP (The GNU Multiple Precision Arithmetic Library [12]) zum Einsatz. Das Programm hat dann auf einem Surface 6 Pro binnen sieben Sekunden die korrekten Lösungen gefunden und dabei nicht eine einzige Falschlösung ausgespuckt. Hierbei waren insgesamt nur acht Threads beteiligt.

int main(int argc, char ** argv) {
    uint64_t const gbegin = std::stoll(argc >= 2 ? argv[1] : "0"),
        gend = std::stoll(argc >= 3 ? argv[2] : "0");
    size_t const nthreads = std::thread::hardware_concurrency();
    std::cout << "Num threads: " << nthreads << std::endl
        << "Begin: " << gbegin << ", End: " << gend << std::endl;

    auto stime = Time();
    if (gbegin < gend) {
        std::mutex cout_mux;
        std::vector<std::future<void>> threads;
        for (size_t ithread = 0; ithread < nthreads; ++ithread)
            threads.emplace_back(std::async(std::launch::async, [&, ithread]{
                uint64_t const block = (gend - gbegin + nthreads - 1) / nthreads,
                    begin = gbegin + ithread * block,
                    end = std::min<size_t>(gbegin + (ithread + 1) * block, gend);
                mpz_class x, t0, t1, y, root, rem;
                for (uint64_t ix = begin; ix < end; ++ix) {
                    x = uint32_t(ix >> 32); x <<= 32; x |= uint32_t(ix);
                    mpz_pow_ui(t0.get_mpz_t(), x.get_mpz_t(), 6);
                    mpz_pow_ui(t1.get_mpz_t(), x.get_mpz_t(), 2);
                    y = t0 - 4 * t1 + 4;
                    mpz_sqrtrem(root.get_mpz_t(), rem.get_mpz_t(), y.get_mpz_t());
                    if (rem == 0) {
                        std::lock_guard<std::mutex> lock(cout_mux);
                        std::cout << "[" << mpz_to_str(x) << ", " << mpz_to_str(root)
                            << ", " << mpz_to_str(y) << "]" << std::endl;
                    }
                }
            }));
    }
    stime = Time() - stime;
    std::cout << "Time " << std::fixed << std::setprecision(3)
        << stime << " sec" << std::endl;
}

Listing 9: Die Routine aus Listing 5 wurde hier in C++ umgesetzt.

JIT, CUDA oder C++ für mehr Performance

Das konsistente rechenintensive (aber nicht komplizierte) Beispiel legt dar, welche Präzisionseinbußen eine wirklich starke Beschleunigung nach sich zieht. Aber auch, mit welchen Performance-Einbußen Entwickler rechnen müssen, wenn Präzision an erster Stelle steht. Zusammenfassend lässt sich sagen, dass Python mittels JIT- oder CUDA-Optimierung zwar eine hoch performante Implementierung erlaubt, Programmierer sich aber in sehr hohen Zahlenregionen mit Überläufen konfrontiert sehen. Der beste Kompromiss zwischen Performance und Präzision ist der Einsatz von gmpy2 (ohne JIT-Optimierung aber mit multiprocessing). Denn in dieser Konstellation sind keine Überläufe und damit einhergehend Falschlösungen entstanden und die Laufzeit betrug etwa 15 Minuten (mit multiprocessing 49 Sekunden) und keine Stunden wie im Falle der (präzisionsverlustfreien) Variante NumPy oder Mathematica.

Als unschlagbar hat sich die Implementierung in C++ herausgestellt: Sie lieferte auf einem Surface 6 Pro nur korrekte Lösungen ohne Präzisionsverlust. Wenn es also ausschließlich um Präzision und Performance geht, werden Programmierer nicht um den Einsatz C oder C++ herumkommen.

Eldar Sultanow
ist Architekt bei Capgemini. Seine Schwerpunkte sind moderne Softwarearchitekturen, Data Science und Unternehmensarchitekturmanagement.

(fms [13])

URL dieses Artikels:
https://www.heise.de/-7127649

Links in diesem Artikel:
[1] https://www.heise.de/tests/GPU-Server-im-iX-Test-Dells-DSS-8440-fuer-KI-und-ML-Anwendungen-4687829.html
[2] https://aws.amazon.com/de/ec2/instance-types/p3/
[3] https://www.heise.de/news/Nvidia-erweitert-das-Programmiermodell-in-CUDA-11-5-6229169.html
[4] https://www.spektrum.de/wissen/pietro-mengoli-1626-1686-ueber-die-unendlichkeit-hinaus/1737066
[5] https://news.mit.edu/2019/answer-life-universe-and-everything-sum-three-cubes-mathematics-0910
[6] https://github.com/Sultanow/pythagorean
[7] https://math.stackexchange.com/a/3278708/993738
[8] https://github.com/Sultanow/elliptic_curves/tree/master/python
[9] https://math.stackexchange.com/q/4347945/993738
[10] https://github.com/Sultanow/elliptic_curves/blob/master/python/find_integer_solutions.py
[11] https://github.com/Sultanow/elliptic_curves/tree/master/cpp
[12] https://gmplib.org/
[13] mailto:fms@heise.de