John Carmacks ungewöhnlich schnelle inverse Quadratwurzel (Quake III)

Question

John Carmack hat eine spezielle Funktion im Quellcode von Quake III, die die inverse Quadratwurzel eines Floats berechnet, 4x schneller als die normale (float)(1.0/sqrt(x)) einschließlich einer seltsamen 0x5f3759df konstant. Siehe den nachstehenden Code. Kann jemand Zeile für Zeile erklären, was genau hier vor sich geht und warum dies so viel schneller funktioniert als die normale Implementierung?

float Q_rsqrt( float number )
{
  long i;
  float x2, y;
  const float threehalfs = 1.5F;

  x2 = number * 0.5F;
  y  = number;
  i  = * ( long * ) &y;
  i  = 0x5f3759df - ( i >> 1 );
  y  = * ( float * ) &i;
  y  = y * ( threehalfs - ( x2 * y * y ) );

  #ifndef Q3_VM
  #ifdef __linux__
    assert( !isnan(y) );
  #endif
  #endif
  return y;
}

Answer 1

ZU IHRER INFORMATION. Carmack hat es nicht geschrieben. Terje Mathisen und Gary Tarolli haben es zum Teil (und in sehr bescheidenem Umfang) geschrieben und auch einige andere Quellen genannt.

Wie die mythische Konstante zustande kam, ist ein Rätsel.

Um Gary Tarolli zu zitieren:

Was eigentlich ein Floating ist Gleitkommaberechnung in Ganzzahl - es dauerte lange Zeit gebraucht, um herauszufinden, wie und warum das funktioniert, und ich kann mich nicht mehr an die Details nicht mehr erinnern.

Eine etwas bessere Konstante, von einem erfahrenen Mathematiker entwickelt (Chris Lomont) versucht, herauszufinden, wie der ursprüngliche Algorithmus funktioniert hat:

float InvSqrt(float x)
{
    float xhalf = 0.5f * x;
    int i = *(int*)&x;              // get bits for floating value
    i = 0x5f375a86 - (i >> 1);      // gives initial guess y0
    x = *(float*)&i;                // convert bits back to float
    x = x * (1.5f - xhalf * x * x); // Newton step, repeating increases accuracy
    return x;
}

Trotzdem erwies sich sein anfänglicher Versuch einer mathematisch "überlegenen" Version von id's sqrt (die fast dieselbe Konstante ergab) als schlechter als die ursprünglich von Gary entwickelte, obwohl sie mathematisch viel "reiner" war. Er konnte nicht erklären, warum id's sqrt so hervorragend war (iirc).

Answer 2

Natürlich stellt sich heutzutage heraus, dass dies viel langsamer ist als die Verwendung von sqrt einer FPU (insbesondere auf 360/PS3), da das Umschalten zwischen Float- und Int-Registern einen Load-Hit-Store verursacht, während die Floating Point Unit reziproke Quadratwurzeln in Hardware ausführen kann.

Es zeigt nur, wie sich Optimierungen weiterentwickeln müssen, wenn sich die Art der zugrunde liegenden Hardware ändert.

Answer 3

Greg Hewgill y IllidanS4 hat einen Link mit einer ausgezeichneten mathematischen Erklärung angegeben. Ich werde versuchen, es hier für diejenigen zusammenzufassen, die nicht zu sehr ins Detail gehen wollen.

Jede mathematische Funktion, mit einigen Ausnahmen, kann durch eine Polynomsumme dargestellt werden:

y = f(x)

kann sein genau in umgewandelt:

y = a0 + a1*x + a2*(x^2) + a3*(x^3) + a4*(x^4) + ...

Dabei sind a0, a1, a2,... Konstanten . Das Problem ist, dass für viele Funktionen, wie z. B. Quadratwurzel, für den exakten Wert diese Summe eine unendliche Anzahl von Gliedern hat, sie endet nicht bei einem x^n . Aber, wenn wir bei einigen x^n hätten wir immer noch ein Ergebnis bis zu einer gewissen Genauigkeit.

Also, wenn wir haben:

y = 1/sqrt(x)

In diesem speziellen Fall haben sie beschlossen, alle Polynomglieder über der zweiten zu verwerfen, wahrscheinlich aus Gründen der Rechengeschwindigkeit:

y = a0 + a1*x + [...discarded...]

Die Aufgabe besteht nun darin, a0 und a1 so zu berechnen, dass y die geringste Abweichung vom exakten Wert aufweist. Sie haben berechnet, dass die am besten geeigneten Werte sind:

a0 = 0x5f375a86
a1 = -0.5

Wenn man dies in die Gleichung einbezieht, erhält man also:

y = 0x5f375a86 - 0.5*x

Das ist die gleiche Zeile wie die im Code:

i = 0x5f375a86 - (i >> 1);

Edit: eigentlich hier y = 0x5f375a86 - 0.5*x ist nicht dasselbe wie i = 0x5f375a86 - (i >> 1); da das Verschieben von Float als Integer nicht nur durch zwei dividiert, sondern auch den Exponenten durch zwei teilt und einige andere Artefakte verursacht, aber es geht immer noch um die Berechnung einiger Koeffizienten a0, a1, a2... .

Zu diesem Zeitpunkt haben sie festgestellt, dass die Genauigkeit dieses Ergebnisses nicht ausreicht. Daher wurde zusätzlich nur ein Schritt der Newtonschen Iteration durchgeführt, um die Genauigkeit des Ergebnisses zu verbessern:

x = x * (1.5f - xhalf * x * x)

Sie hätten einige weitere Iterationen in einer Schleife durchführen können, wobei jede einzelne das Ergebnis verbessert hätte, bis die erforderliche Genauigkeit erreicht ist. Genau so funktioniert es in CPU/FPU! Aber es scheint, dass nur eine Iteration ausreichend war, was auch ein Segen für die Geschwindigkeit war. Die CPU/FPU führt so viele Iterationen durch, wie erforderlich sind, um die Genauigkeit der Fließkommazahl zu erreichen, in der das Ergebnis gespeichert wird, und sie verfügt über einen allgemeineren Algorithmus, der für alle Fälle geeignet ist.

---

Kurz gesagt, was sie getan haben, ist:

Verwenden Sie (fast) denselben Algorithmus wie CPU/FPU, nutzen Sie die Verbesserung der Anfangsbedingungen für den Sonderfall 1/sqrt(x) und rechnen Sie nicht bis zur Genauigkeit, die CPU/FPU erreicht, sondern brechen Sie früher ab, wodurch Sie an Rechengeschwindigkeit gewinnen.

Answer 4

Ich war neugierig zu sehen, was die Konstante als Float war, so schrieb ich einfach dieses Stück Code und googelte die ganze Zahl, die herauskam.

long i = 0x5F3759DF;
float* fp = (float*)&i;
printf("(2^127)^(1/2) = %f\n", *fp);
//Output
//(2^127)^(1/2) = 13211836172961054720.000000

Es sieht so aus, als sei die Konstante "Eine ganzzahlige Annäherung an die Quadratwurzel von 2^127, besser bekannt durch die hexadezimale Form ihrer Gleitkommadarstellung, 0x5f3759df". https://mrob.com/pub/math/numbers-18.html

Auf derselben Website wird das Ganze erklärt. https://mrob.com/pub/math/numbers-16.html#le009_16

Answer 5

Nach zu diesem schönen Artikel die vor einiger Zeit geschrieben wurde...

Die Magie des Codes, auch wenn Sie nicht nachvollziehen kann, sticht die Zeile i = 0x5f3759df - (i>>1); Zeile. Vereinfacht, ist Newton-Raphson eine Näherung die mit einer Schätzung beginnt und mit Iteration verfeinert. Unter Verwendung von Vorteil der Natur von 32-Bit-x86 Prozessoren wird i, eine ganze Zahl, zunächst anfänglich auf den Wert der Gleitkommazahl gesetzt, die Sie inversen Quadrats der Fließkommazahl gesetzt, wobei ein Integer-Cast. i wird dann gesetzt auf 0x5f3759df, minus sich selbst, um ein Bit verschoben Bit nach rechts verschoben. Die Rechtsverschiebung fällt das niedrigstwertige Bit von i, und halbiert es im Wesentlichen.

Es ist eine wirklich gute Lektüre. Dies ist nur ein kleiner Teil davon.