- Die durch Quake 3 bekannt gewordene fast inverse square root war damals eine performante Lösung, um
1 / sqrt(x) durch Neuinterpretation der Float-Bits und eine Newton-Raphson-Korrektur schnell zu approximieren
- Der Kern ist, dass das Integer-Bitmuster eines IEEE-754-32-Bit-Floats wie eine skalierte und verschobene Näherung von
log2(x) behandelt werden kann
0x5f3759df - (i >> 1) ist die auf Integer-Shift und Subtraktion übertragene Form von log2(x^-0.5) = -0.5 * log2(x); die Magic Constant ergibt sich aus 3/2 * 2^23 * (127 - σ)- Anschließend wird mit einmal
y = y * (1.5 - 0.5x * y * y) eine Newton-Raphson-Korrektur angewendet; die zweite Iteration im Quake-Code ist auskommentiert
- 1999 wurde die inverse Quadratwurzel für Beleuchtung und Normalisierung von 3D-Vektoren hunderte bis tausende Male pro Sekunde benötigt; auf moderner Hardware ist der praktische Nutzen desselben Tricks dank dedizierter Gleitkommaeinheiten geringer
Was der Quake-Code macht
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 ) );
return y;
}
- Diese Funktion berechnet einen Näherungswert für die inverse Quadratwurzel
1 / sqrt(number) von number
- Der berühmteste Teil ist die Bitmanipulation, bei der ein Float-Wert wie ein
long interpretiert und dann 0x5f3759df - (i >> 1) ausgeführt wird
- Als Quake 3 im Jahr 1999 erschien, war die inverse Quadratwurzel eine langsame und teure Operation und wurde bei Beleuchtungsgleichungen sowie bei Berechnungen von 3D-Vektoren, die normalisiert werden mussten, hunderte bis tausende Male pro Sekunde benötigt
- Auf moderner Hardware werden solche Berechnungen entweder nicht auf der CPU ausgeführt oder sind selbst auf der CPU dank weiterentwickelter dedizierter Gleitkomma-Hardware schnell
IEEE-754-Darstellung eines 32-Bit-Floats
- Ein 32-Bit-Float besteht aus drei Teilen
- Sign: 1 Bit, gibt an, ob der Wert positiv oder negativ ist
- Exponent: 8 Bit, legt den Wertebereich fest, in dem der Wert liegt
- Mantissa: 23 Bit, gibt linear die Position innerhalb dieses Bereichs an
- Ein normaler Wert wird in folgender Form interpretiert
N = (-1)^S * 2^(E - 127) * (1 + M / 2^23)
B = 127 ist der Bias-Wert für den biased exponent, der tatsächliche Exponent ist e = E - B- Die Mantisse wird nicht einfach mit
m multipliziert, sondern in der Form 1 + m verwendet
- Wenn
m = 0, ergibt sich 2^e
- Wenn
m nahe an 1 heranrückt, werden Werte bis kurz vor den nächsten Exponentenbereich 2^(e+1) dargestellt
- Wenn alle Exponentenbits 0 sind, handelt es sich um eine subnormale Zahl, und die Formel ändert sich
N = (-1)^S * 2^-126 * m
- Subnormale Zahlen sind nötig, um 0 und sehr kleine Zahlen nahe 0 darzustellen
- Wenn alle Exponentenbits 1 sind, wird der Wert als Spezialwert behandelt
- Bei
E = 255, M = 0 ist der Wert Infinity oder -Infinity
- Bei
M != 0 ist der Wert NaN
Die logarithmische Beziehung beim Betrachten von Float-Bits als Integer
- Betrachtet man die interne Darstellung eines Floats wie einen 32-Bit-Integer, lässt sie sich mit folgender Formel ausdrücken
I_x = 2^31 S + 2^23 E + M
- Da die inverse Quadratwurzel positive Eingaben voraussetzt, wird die Formel mit
S = 0 einfacher
L = 2^23
I_x = L E + M
- Innerhalb desselben Exponentenbereichs gibt die Mantisse die Position linear an; mit größerem Exponenten deckt dieselbe Anzahl an Mantissen-Schritten jedoch einen größeren Abschnitt auf der Zahlengeraden ab
E = 127, also e = 0, entspricht ungefähr dem Bereich [1, 2)E = 128, also e = 1, entspricht ungefähr dem Bereich [2, 4)- Beide Bereiche haben dieselbe Anzahl an Mantissen-Schritten, aber der zweite Bereich ist doppelt so breit
- Wegen dieser Struktur entsteht eine logarithmische Beziehung, wenn man das rohe Bitmuster eines Floats als Integer betrachtet
Rohbits als Näherung von log2(x)
- Interpretiert man das Bitmuster eines Floats als Integer
I_x, kann man es als stückweise lineare Näherung von log2(x) betrachten
- Diese Beziehung lässt sich durch folgende Näherungsformel ausdrücken
log2(x) ≈ I_x / L - B
- Teilt man den Rohbit-Integer durch die Mantissengröße
L = 2^23 und zieht den Exponenten-Bias B = 127 ab, erhält man einen Wert nahe log2(x)
- Der Logarithmus innerhalb des Mantissenintervalls wird linear angenähert
log2(1 + x) ≈ x + σ
σ ist ein Tuning-Parameter, der die Näherung justiert, und x gibt die Position innerhalb des Exponentenintervalls im Bereich [0, 1] an
Die inverse Quadratwurzel in eine Logarithmus-Identität umformen
- Ziel ist es, folgenden Wert zu berechnen
y = 1 / sqrt(x)
- In Exponentialform wird daraus
y = x^-0.5
- Wendet man die Logarithmus-Identität an, ergibt sich für die Berechnung der inversen Quadratwurzel folgende Beziehung
log2(1 / sqrt(x)) = log2(x^-0.5) = -0.5 * log2(x)
- Nutzt man aus, dass Float-Bits wie eine Näherung von
log2(x) funktionieren, lässt sich aus der Integer-Bitdarstellung I_x von x die Integer-Bitdarstellung I_y von y direkt approximieren
I_y ≈ -0.5 I_x + 1.5 L (B - σ)
- Diese Formel führt zur zentralen Zeile im Quake-Code
i = 0x5f3759df - ( i >> 1 );
i >> 1 verschiebt die Integer-Bits um 1 Bit nach rechts und wirkt damit wie eine Multiplikation mit 1/2- Die vorangestellte Konstante
0x5f3759df entspricht 1.5 * L * (B - σ)
Was hinter der Konstante 0x5f3759df steckt
- Setzt man
σ = 0, wird die Konstante wie folgt berechnet
1.5 * 2^23 * 127 = 1598029824
- Die hexadezimale Darstellung dieses Werts ist
0x5f400000
- Vom tatsächlichen Quake-Wert
0x5f3759df unterscheidet er sich um 566817
- Aus dieser Differenz lässt sich der zum Quake-Code passende Wert für
σ wie folgt berechnen
σ = 377878 / 2^23
σ = 0.04504656
- In C lässt sich dieselbe Konstante so berechnen
int32_t compute_magic(void) {
double sigma = 0.0450465;
double expression = 1.5 * pow(2.0, 23.0) * (127.0 - sigma);
int32_t i = expression;
return i;
}
// -> 0x5f3759df
- Hier wird
double verwendet, und die Integer-Konvertierung ist keine Bit-Neuinterpretation, sondern ein normaler Cast
- Dieser Wert für
σ wurde gewählt, um die Näherung zu optimieren; er ist jedoch nicht der tatsächlich optimale Wert, und es ist auch nicht sicher, wer ihn erstellt hat
Warum es kein einfacher Hack ist
0x5f3759df - (i >> 1) ist eine Formel, die ausnutzt, dass die Rohbits eines Floats eine Logarithmus-Näherung bilden, um einen Startwert für die inverse Quadratwurzel zu erzeugen- Sie basiert auf komplexen mathematischen Beziehungen, verwendet zur Laufzeit aber nur schnelle Operationen wie Shift und Subtraktion
- Da damals teure Operationen tausende Male pro Sekunde ausgeführt werden mussten, wurde diese Methode zu einem an die Hardwarebeschränkungen angepassten Engineering-Design
- Allerdings funktioniert dieser Algorithmus nur mit normalen Floats
- Bei subnormalen Werten gilt die Annahme der Näherung
log2(1 + x) ≈ x + σ nicht
- Bei subnormalen Werten steckt tatsächlich eher eine Form nahe
0 + x darin, wodurch die Näherung zusammenbricht
Fehlerreduktion durch Newton-Raphson-Korrektur
- Der durch die Bitmanipulation gewonnene Startwert ist ziemlich gut, enthält aber noch messbare Fehler
- Die folgende Zeile verbessert die Näherung deutlich
y = y * ( threehalfs - ( x2 * y * y ) );
- Diese Zeile ist eine Anwendung der Newton-Raphson method
- Um das Problem der inversen Quadratwurzel an die Newton-Methode anzupassen, wird es in die Suche nach einer Nullstelle der folgenden Funktion umgeformt
f(y) = 1 / y^2 - x = 0
- Die Newton-Methode erzeugt aus dem aktuellen Näherungswert
y_n wie folgt einen besseren Näherungswert y_(n+1)
y_(n+1) = y_n - f(y_n) / f'(y_n)
- Die Ableitung von
f(y) = y^-2 - x lautet
f'(y) = -2y^-3 = -2 / y^3
Newton-Korrekturformel ohne Division
- Wendet man die Newton-Formel direkt an, enthält sie mehrere Gleitkomma-Divisionen
- Einer der Gründe für die Geschwindigkeit dieses Algorithmus ist die Vermeidung von Gleitkomma-Divisionen
- Algebraisch umgeformt entsteht eine Variante, die ohne Division auskommt und nur Multiplikationen nutzt
y_(n+1) = y_n * (1.5 - 0.5x * y_n^2)
- Im Quake-Code wird mit
x2 = number * 0.5F der Wert 0.5x vorab berechnet und in der folgenden Zeile verwendet
y = y * ( threehalfs - ( x2 * y * y ) );
- Nach dieser einen Iteration beträgt der maximale absolute Fehler 0,175 %, und in vielen Fällen ist der Fehler kleiner
- Im Originalcode gibt es eine zweite Newton-Iteration, sie ist jedoch auskommentiert
// y = y * ( threehalfs - ( x2 * y * y ) ); // 2nd iteration, this can be removed
Ursprung und verwandte Algorithmen
- Dieser Algorithmus wurde nicht von John Carmack erfunden, und sein genauer Ursprung ist nicht zu 100 % sicher
- Dazu ist ein Artikel von Beyond3D verlinkt: The truth is the exact origin is not 100% certain
- Chris Lomont verfasste ein Paper, das im Schritt der Logarithmus-Näherung den optimalen Sigma-Wert sucht: InvSqrt.pdf
- CORDIC ist ein Algorithmus, der Sinus und Kosinus ohne Gleitkomma, nur mit Additionen und Bit-Shifts berechnet; in den Details unterscheidet er sich stark von fast inverse square root
- Gemeinsam ist beiden Algorithmen, dass sie mathematische Beobachtungen effizient auf die damaligen Hardwarebeschränkungen anwenden
2 Kommentare
Dieser faszinierende Code taucht immer wieder auf, sobald man ihn fast vergessen hat.. haha
Meinungen auf Hacker News
Wenn ein Computer nach 1999 gebaut wurde, unterstützt er in der Regel den SSE-Befehlssatz; dazu gehört
_mm_rsqrt_ps, das vier inverse Quadratwurzeln auf einmal schneller berechnet: https://www.intel.com/content/www/us/en/docs/intrinsics-guid...Trotzdem ist die hier behandelte Technik noch nicht völlig bedeutungslos. Float/Int-Konvertierungen sind schnell, aber es gibt weiterhin Hardware ohne
rsqrt-,sqrt-,pow- oderlog-Befehle, und solche Operationen lassen sich mit diesem Trick approximieren.sqrt(x)alsx * 1/sqrt(x)zu berechnen.In GPU-Befehlssätzen, ARM, RISC-V, AVR, PIC, 8051, FPGAs usw. sind approximative inverse Quadratwurzel-Operationen oft eingebaut; vermutlich sind sie aber mit solchen Algorithmen implementiert.
Wenn man am Artikel etwas herumkritteln will: Die Erklärung, solche Berechnungen würden auf heutigen CPUs nicht mehr stattfinden, stimmt so nicht. Es ist ein verbreiteter Irrtum, dass Spiele oder Apps mit vielen Gleitkommaoperationen alle Gleitkommaoperationen an die GPU auslagern wollen.
Sinnvoll ist das Auslagern an die GPU in der Praxis nur bei großen, gleichförmigen Aufgaben. Wenn man eine einmalige Vektornormalisierung macht, etwa um eine Rotationsmatrix aufzubauen, damit ein Objekt ein anderes anschaut, ist es schneller, das auf der CPU zu lassen. Selbst ohne die Übertragungszeit zur GPU ist eine einzelne Gleitkommaoperation auf der CPU schneller, weil GPUs in der Regel niedriger takten und ihre hohen FLOP-Zahlen durch Parallelität erreichen.
Ich habe eine MMIX-Implementierung geschrieben und setze dabei voraus, dass die ursprüngliche Eingabe größer als
2^-1021ist.Falls es interessiert: Auch Wikipedia hat eine recht gute Erklärung dieser Funktion und ihrer Geschichte: https://en.wikipedia.org/wiki/Fast_inverse_square_root
Ich habe ein paar solcher Dinge gesammelt: https://github.com/ncruces/fastmath/blob/main/fast.go
Dazu gibt es auch einen passenden StackOverflow-Beitrag: https://stackoverflow.com/questions/32042673/optimized-low-a...
fastmath-Pakets sehen.Zeit zum Korinthenkacken. In der Float-Formel gibt es einen Tippfehler: Es muss
(-1)^Sheißen, nicht-1^S. Letzteres ergibt immer-1.Auch die Erklärung, das Interpretieren des rohen Bitmusters sei eine stückweise lineare Approximation des Logarithmus, ist nicht ganz korrekt. Die Linien zwischen den Datenpunkten im blauen Graphen existieren tatsächlich nicht, und ein Bit kann nicht nur zur Hälfte auf 1 gesetzt sein. Es ähnelt eher einer diskreten Version des Logarithmus, und die tatsächlich existierenden Datenpunkte – also die Punkte, an denen die rote und die blaue Linie zusammentreffen – sind buchstäblich gleich einem skalierten und verschobenen Logarithmus. Ansonsten ein guter Text.
[010000, 010111]die Werte 2, 2,25, 2,5, 2,75, 3, 3,25, 3,5, 3,75.Die durch den Logarithmus zur Basis 2 dieser Zahlen implizierten Mantissen sind aber jeweils
.0000000,.0010101,.0101001,.0111010,.1001010,.1011001,.1100111,.1110100; abgesehen vom ersten Wert entsprechen sie nicht den Float-Mantissen001,010usw. Da Floats im Intervall[2,4)linear beabstandet sind, die entsprechenden Logarithmen aber nicht, kann man Floats, wie der Artikel sagt, als stückweise lineare Approximation des Logarithmus betrachten.Für den vollständigen Graphen gäbe es innerhalb des stückweise linearen Musters
2^32Auswahlmöglichkeiten, aber das, was der Originalbeitrag zeichnet, ist nicht dieser vollständige Graph. Da der Text 32-Bit-Integer und IEEE-754-32-Bit-Float-Operationen behandelt, finde ich es in Ordnung, in der Erklärung „diskret“ wegzulassen.Ein guter Text, der viele interessante Konzepte erklärt, aber die algebraische Herleitung in einem Abschnitt ist erstaunlich schlecht.
Nach „Die genauen Schritte, um von der ersten Form zu dieser Form zu gelangen, sind zahlreich, aber der Vollständigkeit halber sind sie alle enthalten“ enthält die Herleitung viele unnötige Schritte und mehrere Vorzeichenfehler, die sich gegenseitig aufheben. Insbesondere wird beim Übergang von der zweiten zur dritten Zeile das Minuszeichen nicht korrekt verteilt. Nach der zweiten Zeile könnte man bei
y_n+1 = y_n + (1 - x * y_n^2) / y_n^2 * (y_n^3 / 2)beginnen und viel kürzer zuy_n+1 = y_n (1.5 * y_n - 0.5 * x * y_n * y_n)gelangen; auch die Zwischenschritte wären dann korrekt. Für jemanden, der Algebra versteht, sind das meiner Meinung nach nur offensichtliche Schritte.Die Magic Number des berühmten Code-Schnipsels ist nicht die optimale Konstante. Mit einer anderen Konstante lässt sich der relative Fehler vermutlich um etwa 0,5 % weiter verringern.
Damals war es vielleicht schwer, den absolut optimalen Wert zu finden, heute ist es vergleichsweise einfach. Ich bin selbst einmal in dieses Rabbit Hole geraten und habe ein Jupyter-Notebook, das die optimalen Magic Numbers für
(1/x^2)und(1/x)sucht.Das Interessanteste an diesem Artikel war für mich der Link zu „How Java's Floating-Point Hurts Everyone Everywhere“: https://people.eecs.berkeley.edu/~wkahan/JAVAhurt.pdf
Der Autor ist William Kahan, auch bekannt als „Old Man of Floating-Point“: https://news.ycombinator.com/item?id=29042853 – An Interview with the Old Man of Floating-Point (1998)
JAVAhurt-PDF zu lesen, und der Satz ist schrecklich. Es wirkt, als wäre ein TeX-Paket verwendet worden, das Wortabstände übermäßig und dazu noch ungleichmäßig auseinanderzieht, oder als wäre ein anderes Dokument per OCR erfasst worden und hätte zusätzliche Leerzeichen bekommen.Auch in den Abschnitten mit Monospace-Schrift gibt es merkwürdige zusätzliche Abstände. Es war wirklich schwer, sich beim Lesen zu konzentrieren, und obwohl ich weiß, dass es das nicht ist, fühlte es sich fast wie ein Manifest eines Wissenschafts-Cranks an.
Dieses Video, das ich vor einiger Zeit gesehen habe, fand ich wirklich gut: https://www.youtube.com/watch?v=p8u_k2LIZyo