John Carmacks ungewöhnliche schnelle Quadratwurzel (Quake III)

112

John Carmack hat eine spezielle Funktion im Quake III-Quellcode, die die inverse Quadratwurzel eines Floats berechnet, 4x schneller als normal (float)(1.0/sqrt(x)), einschließlich einer seltsamen 0x5f3759dfKonstante. Siehe den Code unten. Kann jemand Zeile für Zeile erklären, was genau hier vor sich geht und warum dies so viel schneller funktioniert als die reguläre 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;
}
Alex
quelle
10
Dies wurde über zig Mal geschrieben. Siehe: google.com/search?q=0x5f3759df
Greg Hewgill
15
Trotzdem danke. Dies war eine viel interessantere Frage als "Wie macht man eine positive Zahl in C # negativ?"
MusiGenesis
7
Heiliger Mist, dies ist nur ein Hack, der auf Newtons Methode basiert, es ist kein heiliger Gral von Algorithmen, hör auf darüber zu reden,
bitte

Antworten:

75

Zu Ihrer Information. Carmack hat es nicht geschrieben. Terje Mathisen und Gary Tarolli nehmen beide teilweise (und sehr bescheiden) Anerkennung dafür sowie einige andere Quellen.

Wie die mythische Konstante abgeleitet wurde, ist ein Rätsel.

Um Gary Tarolli zu zitieren:

Was eigentlich eine Gleitkomma-Berechnung in Ganzzahl macht - es hat lange gedauert, um herauszufinden, wie und warum dies funktioniert, und ich kann mich nicht mehr an die Details erinnern.

Eine etwas bessere Konstante, die von einem erfahrenen Mathematiker (Chris Lomont) entwickelt wurde, um herauszufinden, wie der ursprüngliche Algorithmus funktioniert, ist:

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 erster Versuch, eine mathematisch "überlegene" Version von ids sqrt (die fast dieselbe Konstante erreichte), als schlechter als die ursprünglich von Gary entwickelte, obwohl er mathematisch viel "reiner" war. Er konnte nicht erklären, warum es so gut war.

Rushyo
quelle
4
Was soll "mathematisch reiner" bedeuten?
Tara
1
Ich würde mir vorstellen, wo die erste Vermutung aus berechtigten Konstanten abgeleitet werden kann, anstatt scheinbar willkürlich zu sein. Wenn Sie eine technische Beschreibung wünschen, können Sie diese nachschlagen. Ich bin kein Mathematiker und eine semantische Diskussion über mathematische Terminologie gehört nicht zu SO.
Rushyo
7
Das ist genau der Grund, warum ich dieses Wort in Angstzitaten zusammengefasst habe, um diese Art von Unsinn abzuwenden. Das setzt voraus, dass der Leser mit der umgangssprachlichen englischen Schrift vertraut ist, denke ich. Sie würden denken, dass gesunder Menschenverstand ausreichen würde. Ich habe keinen vagen Begriff verwendet, weil ich dachte: "Weißt du was? Ich möchte wirklich von jemandem befragt werden, der sich nicht die Mühe macht, die Originalquelle nachzuschlagen, was bei Google zwei Sekunden dauern würde."
Rushyo
2
Nun, Sie haben die Frage tatsächlich nicht beantwortet.
BJovke
1
Für diejenigen, die wissen wollten, wo er es findet: beyond3d.com/content/articles/8
mr5
52

Natürlich stellt sich heutzutage heraus, dass es viel langsamer ist als nur die Verwendung des sqrt einer FPU (insbesondere unter 360 / PS3), da das Wechseln zwischen float- und int-Registern einen Load-Hit-Store induziert, während die Gleitkommaeinheit ein reziprokes Quadrat ausführen kann root in Hardware.

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

Crashworks
quelle
4
Es ist immer noch viel schneller als std :: sqrt ().
Tara
2
Hast du eine Quelle? Ich möchte die Laufzeiten testen, habe aber kein Xbox 360-Entwicklungskit.
DucRP
31

Greg Hewgill und IllidanS4 gaben einen Link mit ausgezeichneter mathematischer Erklärung. Ich werde versuchen, es hier für diejenigen zusammenzufassen, die nicht zu sehr ins Detail gehen wollen.

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

y = f(x)

kann genau umgewandelt werden in:

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

Wobei a0, a1, a2, ... Konstanten sind . Das Problem ist, dass für viele Funktionen, wie die Quadratwurzel, diese Summe für den exakten Wert unendlich viele Mitglieder hat und nicht bei x ^ n endet . Aber wenn wir bei x ^ n anhalten haben wir immer noch ein Ergebnis mit einer gewissen Präzision.

Also, wenn wir haben:

y = 1/sqrt(x)

In diesem speziellen Fall haben sie beschlossen, alle Polynomelemente über der Sekunde zu verwerfen, wahrscheinlich aufgrund der Berechnungsgeschwindigkeit:

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

Und jetzt ist die Aufgabe gekommen, a0 und a1 zu berechnen, damit y den geringsten Unterschied zum exakten Wert aufweist. Sie haben berechnet, dass die am besten geeigneten Werte sind:

a0 = 0x5f375a86
a1 = -0.5

Wenn Sie dies in eine Gleichung setzen, erhalten Sie:

y = 0x5f375a86 - 0.5*x

Welches ist das gleiche wie die Zeile, die Sie im Code sehen:

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

Edit: eigentlich ist hier y = 0x5f375a86 - 0.5*xnicht das gleiche wiei = 0x5f375a86 - (i >> 1); seit dem Verschieben von Float als Ganzzahl nicht nur durch zwei geteilt, sondern auch Exponent durch zwei geteilt und einige andere Artefakte verursacht, aber es kommt immer noch darauf an, einige Koeffizienten a0, a1, a2 ... zu berechnen.

Zu diesem Zeitpunkt haben sie herausgefunden, dass die Genauigkeit dieses Ergebnisses für diesen Zweck nicht ausreicht. Daher haben sie zusätzlich nur einen Schritt der Newtonschen Iteration durchgeführt, um die Ergebnisgenauigkeit zu verbessern:

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

Sie hätten einige weitere Iterationen in einer Schleife durchführen können, wobei jede das Ergebnis verbessert, bis die erforderliche Genauigkeit erreicht ist. Genau so funktioniert es in CPU / FPU! Aber es scheint, dass nur eine Iteration ausreichte, was auch ein Segen für die Geschwindigkeit war. Die CPU / FPU führt so viele Iterationen wie nötig durch, um die Genauigkeit für die Gleitkommazahl zu erreichen, in der das Ergebnis gespeichert ist, und verfügt über einen allgemeineren Algorithmus, der in allen Fällen funktioniert.


Kurz gesagt, was sie getan haben, ist:

Verwenden Sie (fast) den gleichen Algorithmus wie CPU / FPU, nutzen Sie die Verbesserung der Anfangsbedingungen für den Sonderfall 1 / sqrt (x) und berechnen Sie nicht bis zur Präzision, zu der CPU / FPU gehen, sondern früher aufhören Berechnungsgeschwindigkeit gewinnen.

BJovke
quelle
2
Das Umsetzen des Zeigers auf long ist eine Annäherung an log_2 (float). Das Zurückwerfen ist eine Annäherung von 2 ^ lang. Dies bedeutet, dass Sie das Verhältnis ungefähr linear machen können.
wizzwizz4
22

Nach diesem schönen Artikel vor einiger Zeit geschrieben ...

Die Magie des Codes, auch wenn Sie ihm nicht folgen können, sticht hervor als i = 0x5f3759df - (i >> 1); Linie. Vereinfacht ausgedrückt ist Newton-Raphson eine Näherung, die mit einer Vermutung beginnt und diese durch Iteration verfeinert. Unter Ausnutzung der Natur von 32-Bit-x86-Prozessoren wird i, eine Ganzzahl, zunächst mit einer Ganzzahlumwandlung auf den Wert der Gleitkommazahl gesetzt, deren inverses Quadrat Sie verwenden möchten. i wird dann auf 0x5f3759df gesetzt, minus selbst um ein Bit nach rechts verschoben. Die Rechtsverschiebung lässt das niedrigstwertige Bit von i fallen und halbiert es im Wesentlichen.

Es ist eine wirklich gute Lektüre. Dies ist nur ein winziges Stück davon.

Dillie-O
quelle
19

Ich war neugierig zu sehen, was die Konstante als Float war, also schrieb ich einfach dieses Stück Code und googelte die Ganzzahl, die heraussprang.

    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 ob die Konstante "Eine ganzzahlige Annäherung an die Quadratwurzel von 2 ^ 127 ist, besser bekannt durch die hexadezimale Form ihrer Gleitkommadarstellung 0x5f3759df" https://mrob.com/pub/math/numbers-18.html

Auf der gleichen Seite erklärt es das Ganze. https://mrob.com/pub/math/numbers-16.html#le009_16

ThisIsAReallyOldQuestion
quelle
6
Dies verdient mehr Aufmerksamkeit. Es macht alles Sinn, nachdem man erkannt hat, dass es nur die Quadratwurzel von 2 ^ 127 ist ...
u8y7541