Zahlen, bitte! 0x5f3759df – merkwürdige Mathematik im Egoshooter

Passend zur Veröffentlichung der neuesten Ausgabe des "Auf dem Mars bricht die Hölle los"-Shooters Doom geht es diesmal um einen Programmiertrick zur besonders schnellen Berechnung von 1/sqrt(x), bei dem die Zahl 0x5f3759df eine unerwartete Rolle spielt. Viele rechnen ihn weiterhin dem Doom-Schöpfer John Carmack zu, weil er im Quellcode von "Quake III Arena" steckte. Tatsächlich ist der Kniff jedoch älter und stammt vom Computergrafik-Veteranen Greg Walsh.

Schnelle Berechnung von Flächennormalen

Zahlen, bitte!

In dieser Rubrik stellen wir immer dienstags verblüffende, beeindruckende, informative und witzige Zahlen aus den Bereichen IT, Wissenschaft, Kunst, Wirtschaft, Politik und natürlich der Mathematik vor.

• Alle Artikel zu "Zahlen, bitte!"

Aber wo braucht man bei Quake solche Berechnungen und warum wurden nicht die Standardfunktionen dafür verwendet? Bei 3D-Grafikprogrammierung benötigt man sie zur Berechnung des Einheitsvektors, der senkrecht auf einer Oberfläche steht (Flächennormale), etwa um Reflexionseigenschaften für Beleuchtungseffekte zu berechnen, statt mit statischen Lightmaps zu arbeiten. Zur Normierung dividiert man den Vektor durch seine Länge – die Wurzel aus der Summe der Quadrate seiner Komponenten ||v|| = sqrt(v_x² + v_y² + v_z²).

Da die Berechnung der Wurzel auf den Prozessoren Ende der 90er-Jahre viel zu langsam gewesen wäre, mussten die Entwickler tief in die Trickkiste greifen. Tatsächlich lief der folgende Code sogar schneller als das sonst übliche Nachschlagen der Ergebnisse in einer Lookup Table:

float Q_rsqrt( float number )
{
long i;
float x2, y;
const float threehalfs = 1.5F;
x2 = number * 0.5F;
y = number;
i = * ( long * ) &y; // evil floating point bit level hacking
i = 0x5f3759df - ( i >> 1 ); // what the fuck?
y = * ( float * ) &i;
y = y * ( threehalfs - ( x2 * y * y ) ); // 1st iteration
// y = y * ( threehalfs - ( x2 * y * y ) ); // 2nd iteration, this can be removed
return y;
}

Auf den ersten Blick käme wohl niemand auf die Idee, dass die Funktion Q_rsqrt() die gesuchte reziproke Quadratwurzel berechnet. Sie wird allerdings auch nicht nicht exakt berechnet, sondern mit Hilfe des Newton-Verfahrens in erster Näherung bestimmt; genauer brauchten es die Entwickler nicht (siehe den Hinweis im Quellcode). Mit dem Verfahren bestimmt man näherungsweise die Nullstellen der Funktion f, wobei f' die Ableitung der Funktion f ist:

yn+1 = yn - f(yn)/f'(yn)

Setzt man f(y) = 1/y² - x, dann findet man die positive Quadratwurzel von x:

Daraus ergibt sich dann:

yn+1 = yn(1.5 − 0.5 * x * yn2)

oder nach Ersetzen der Variablen der oben im Listing gezeigte Code ... so weit, so gut.

Darstellungsfrage

Um schnell zu einer guten Näherung zu kommen, braucht man einen passenden Startwert für die Iteration. Am besten wäre es, dabei zeitaufwendige Operationen wie Divisionen zu vermeiden.

Das "What the fuck?" besteht darin, im Computer nach IEEE 754-1985 dargestellten 32-Bit-Gleitkommazahlen an geeigneter Stelle als Integer zu interpretieren. Diese bestehen aus einem Vorzeichen S (Bit 31), Exponent E (Bit 30 bis 23) und Mantisse M (Bit 22 bis 0) und sehen folgendermaßen aus: y = (-1)^S * (1.M) * 2^E-B, mit B = 127 als Bias (siehe Wikipedia zu Gleitkommazahlen). In Integer-Darstellung ergibt sich i = S * 2³¹ + E * 2²³ + M. Da man für die Berechnung der Flächennormalen nur die positiven Ergebnisse braucht, ist S = 0, was die ganze Sache vereinfacht.

Der Rest besteht in Logarithmieren, Näherungen und der Verwendung der Integer-Repräsentation. Vereinfacht man das Ergebnis, kommt man schließlich zu der Beziehung Iy = R - Ix/2.Das entspricht der Form der Zeile i = 0x5f3759df - ( i >> 1 ), da der Shift-Operator >> die Zahl um ein Bit nach rechts schiebt, was einer Ganzzahldivision durch 2 entspricht, sowohl bei Mantisse als auch beim Exponenten. Danach wird der Wert als Float interpretiert und in das Newton-Verfahren gesteckt.

Die Berechnung des tatsächlichen Werts von R sprengt den Rahmen des Artikels. Wer sich dafür interessiert, findet die Herleitung beispielsweise in "The Mathematics behind the Fast Inverse Square Roote Function Code" von Charles McEniry (PDF).

Heutzutage gibt es für solche Berechnungen spezielle Befehle in CPUs und GPUs, sodass man zumindest auf diesen Trick nicht mehr angewiesen ist. (vza)