- Le fast inverse square root, rendu célèbre par Quake 3, était à l’époque une solution de performance permettant d’approximer rapidement
1 / sqrt(x) au moyen d’une réinterprétation des bits d’un float et d’une correction de Newton-Raphson
- L’idée centrale est que le motif binaire entier d’un float 32 bits IEEE-754 peut être traité comme une approximation mise à l’échelle et décalée de
log2(x)
0x5f3759df - (i >> 1) transpose log2(x^-0.5) = -0.5 * log2(x) en décalage entier et soustraction, et la constante magique vient de 3/2 * 2^23 * (127 - σ)- Ensuite, une seule itération
y = y * (1.5 - 0.5x * y * y) applique la correction de Newton-Raphson ; la seconde itération du code de Quake est commentée
- En 1999, la racine carrée inverse était nécessaire des centaines à des milliers de fois par seconde pour l’éclairage et la normalisation de vecteurs 3D, mais sur le matériel moderne l’intérêt pratique de cette astuce a diminué grâce aux unités de calcul en virgule flottante dédiées
Ce que fait le code de Quake
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;
}
- Cette fonction calcule une approximation de la racine carrée inverse
1 / sqrt(number) pour number
- La partie la plus célèbre est la manipulation de bits qui interprète la valeur float comme un
long, puis exécute 0x5f3759df - (i >> 1)
- En 1999, quand Quake 3 est sorti, la racine carrée inverse était une opération lente et coûteuse, nécessaire des centaines à des milliers de fois par seconde dans les calculs de vecteurs 3D pour les équations d’éclairage et la normalisation
- Sur le matériel moderne, ces calculs ne sont pas effectués par le CPU, ou bien, même lorsqu’ils le sont, ils sont rapides grâce aux unités matérielles dédiées au calcul en virgule flottante
Représentation d’un float 32 bits IEEE-754
- Un float 32 bits est composé de trois parties
- Sign : 1 bit, indiquant si la valeur est positive ou négative
- Exponent : 8 bits, définissant la plage à laquelle appartient la valeur
- Mantissa : 23 bits, représentant linéairement la position dans cette plage
- Une valeur ordinaire est interprétée sous la forme suivante
N = (-1)^S * 2^(E - 127) * (1 + M / 2^23)
B = 127 est la valeur de biais utilisée pour le biased exponent, et l’exposant réel est e = E - B- La mantisse n’est pas simplement utilisée en multipliant par
m, mais sous la forme 1 + m
- si
m = 0, on obtient 2^e
- lorsque
m s’approche de 1, on représente des valeurs juste avant la plage d’exposant suivante, 2^(e+1)
- Si tous les bits de l’exposant valent 0, il s’agit d’un nombre sub-normal, et la formule change
N = (-1)^S * 2^-126 * m
- Les sub-normaux sont nécessaires pour représenter 0 et les très petites valeurs proches de 0
- Si tous les bits de l’exposant valent 1, la valeur est traitée comme une valeur spéciale
- si
E = 255, M = 0, il s’agit de Infinity ou -Infinity
- si
M != 0, il s’agit de NaN
La relation logarithmique qui apparaît quand on voit les bits d’un float comme un entier
- Si l’on voit la représentation interne d’un float comme un entier 32 bits, elle peut être exprimée par la formule suivante
I_x = 2^31 S + 2^23 E + M
- La racine carrée inverse visant des entrées positives, poser
S = 0 simplifie la formule
L = 2^23
I_x = L E + M
- Dans une même plage d’exposant, la mantisse représente linéairement la position, mais plus l’exposant augmente, plus un même nombre de pas de mantisse couvre un intervalle large sur la droite des nombres
E = 127, c’est-à-dire e = 0, correspond approximativement à la plage [1, 2)E = 128, c’est-à-dire e = 1, correspond approximativement à la plage [2, 4)- les deux plages ont le même nombre de pas de mantisse, mais la seconde est deux fois plus large
- À cause de cette structure, voir le motif binaire brut d’un float comme un entier fait apparaître une relation logarithmique
Les bits bruts sont une approximation de log2(x)
- Si l’on interprète le motif binaire d’un float comme un entier
I_x, on peut le voir comme une approximation linéaire par morceaux de log2(x)
- Cette relation s’exprime par l’approximation suivante
log2(x) ≈ I_x / L - B
- Diviser l’entier des bits bruts par la taille de la mantisse
L = 2^23 et soustraire le biais d’exposant B = 127 donne une valeur proche de log2(x)
- Le logarithme à l’intérieur d’une plage de mantisse est traité par approximation linéaire
log2(1 + x) ≈ x + σ
σ est un paramètre de réglage qui ajuste l’approximation, et x indique la position dans la plage d’exposant, dans l’intervalle [0, 1]
Transformer la racine carrée inverse avec une identité logarithmique
- L’objectif est de calculer la valeur suivante
y = 1 / sqrt(x)
- Elle se réécrit sous forme exponentielle comme suit
y = x^-0.5
- En appliquant l’identité logarithmique, le calcul de la racine carrée inverse devient la relation suivante
log2(1 / sqrt(x)) = log2(x^-0.5) = -0.5 * log2(x)
- En exploitant le fait que les bits d’un float se comportent comme une approximation de
log2(x), on peut approximer directement la représentation entière des bits de y, I_y, à partir de la représentation entière des bits de x, I_x
I_y ≈ -0.5 I_x + 1.5 L (B - σ)
- Cette formule mène à la ligne centrale du code de Quake
i = 0x5f3759df - ( i >> 1 );
i >> 1 joue le rôle d’un décalage d’un bit vers la droite sur les bits entiers, ce qui revient à multiplier par 1/2- La constante placée devant,
0x5f3759df, correspond à 1.5 * L * (B - σ)
La véritable nature de la constante 0x5f3759df
- Si l’on pose
σ = 0, la constante se calcule ainsi
1.5 * 2^23 * 127 = 1598029824
- Sa représentation hexadécimale est
0x5f400000
- Elle diffère de la constante réelle de Quake,
0x5f3759df, de 566817
- À partir de cette différence, on peut calculer la valeur de
σ correspondant au code de Quake comme suit
σ = 377878 / 2^23
σ = 0.04504656
- En C, on peut calculer la même constante comme suit
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
- Ici, on utilise
double, et la conversion en entier est un cast ordinaire, pas une réinterprétation des bits
- Cette valeur de
σ a été choisie pour optimiser l’approximation, mais ce n’est pas la vraie valeur optimale, et on ne sait pas non plus avec certitude qui l’a créée
Pourquoi ce n’est pas un simple hack
0x5f3759df - (i >> 1) est une formule qui exploite le fait que les bits bruts d’un float sont une approximation logarithmique afin de produire une valeur initiale pour la racine carrée inverse- Elle repose sur une relation mathématique complexe, mais à l’exécution elle n’utilise que des opérations rapides comme un décalage et une soustraction
- À l’époque, il fallait traiter des opérations coûteuses des milliers de fois par seconde, et cette approche était donc une conception d’ingénierie adaptée aux contraintes matérielles
- Cet algorithme ne fonctionne toutefois que pour les normal float
- pour les valeurs sub-normales, l’hypothèse d’approximation
log2(1 + x) ≈ x + σ ne tient pas
- avec les sub-normaux, on a en réalité une forme proche de
0 + x, ce qui casse l’approximation
Réduire l’erreur avec la correction de Newton-Raphson
- La valeur initiale obtenue par manipulation de bits est assez bonne, mais il reste une erreur mesurable
- La ligne suivante améliore fortement l’approximation
y = y * ( threehalfs - ( x2 * y * y ) );
- Cette ligne applique la Newton-Raphson method
- Pour adapter le problème de la racine carrée inverse à la méthode de Newton, on le transforme en recherche de racine de la fonction suivante
f(y) = 1 / y^2 - x = 0
- La méthode de Newton produit une meilleure approximation
y_(n+1) à partir de l’approximation actuelle y_n, comme suit
y_(n+1) = y_n - f(y_n) / f'(y_n)
- La dérivée de
f(y) = y^-2 - x est la suivante
f'(y) = -2y^-3 = -2 / y^3
Formule de correction de Newton sans division
- Écrire directement la formule de Newton introduit plusieurs divisions en virgule flottante
- L’une des raisons pour lesquelles cet algorithme est rapide est qu’il évite les divisions en virgule flottante
- En la réarrangeant algébriquement, on obtient une forme qui n’utilise que des multiplications, sans division
y_(n+1) = y_n * (1.5 - 0.5x * y_n^2)
- Dans le code de Quake,
x2 = number * 0.5F précalcule 0.5x, qui est ensuite utilisé dans la ligne suivante
y = y * ( threehalfs - ( x2 * y * y ) );
- Après cette unique itération, l’erreur absolue maximale est de 0,175 %, et dans de nombreux cas l’erreur est inférieure à cela
- Le code original contient une seconde itération de Newton, mais elle est commentée
// y = y * ( threehalfs - ( x2 * y * y ) ); // 2nd iteration, this can be removed
Origine et algorithmes liés
- Cet algorithme n’a pas été inventé par John Carmack, et son origine exacte n’est pas connue avec une certitude absolue
- Un article de Beyond3D est lié à ce sujet : The truth is the exact origin is not 100% certain
- Chris Lomont a rédigé un article cherchant la valeur optimale de sigma à l’étape de l’approximation logarithmique : InvSqrt.pdf
- CORDIC est un algorithme qui calcule sine et cosine uniquement avec des additions et des décalages de bits, sans virgule flottante, et son fonctionnement détaillé diffère fortement du fast inverse square root
- Les deux algorithmes ont en commun d’appliquer efficacement des observations mathématiques aux contraintes matérielles de leur époque
2 commentaires
Un bout de code étonnant qui refait surface de temps en temps... haha
Avis sur Hacker News
La plupart des ordinateurs fabriqués depuis 1999 prennent en charge le jeu d’instructions SSE, qui inclut
_mm_rsqrt_pspour calculer plus rapidement quatre inverses de racines carrées à la fois : https://www.intel.com/content/www/us/en/docs/intrinsics-guid...Cela dit, la technique présentée ici n’est pas encore totalement dénuée d’intérêt. Les conversions float/int sont rapides, mais il existe encore du matériel dépourvu d’instructions
rsqrt,sqrt,powoulog, et ces opérations peuvent être approximées avec cette astuce.sqrt(x)commex * 1/sqrt(x).Les jeux d’instructions GPU, ARM, RISC-V, AVR, PIC, 8051, les FPGA, etc. intègrent souvent une opération d’inverse de racine carrée approximative, mais il est probable qu’elle soit implémentée avec ce type d’algorithme.
Pour pinailler un peu sur l’article, l’explication suggérant que ce genre de calcul n’a plus lieu sur les CPU modernes n’est pas correcte. C’est une idée reçue courante de penser que les jeux ou les applications lourdes en calculs en virgule flottante veulent envoyer toutes leurs opérations flottantes au GPU.
En pratique, seuls les gros traitements uniformes ont intérêt à être déportés sur le GPU. Pour une normalisation ponctuelle de vecteur, par exemple construire une matrice de rotation afin qu’un objet en regarde un autre, il est plus rapide de rester sur le CPU. Même sans compter le temps de transfert vers le GPU, une opération flottante isolée est plus rapide sur CPU, car les GPU ont généralement une fréquence plus basse et obtiennent leur nombre élevé de FLOP par le parallélisme.
J’ai écrit une implémentation MMIX, en supposant que l’entrée initiale est supérieure à
2^-1021.Si cela vous intéresse, Wikipedia propose aussi une bonne explication de cette fonction et de son histoire : https://en.wikipedia.org/wiki/Fast_inverse_square_root
J’ai rassemblé quelques éléments de ce genre ici : https://github.com/ncruces/fastmath/blob/main/fast.go
Il existe aussi un fil StackOverflow connexe : https://stackoverflow.com/questions/32042673/optimized-low-a...
fastmath.L’heure est venue de pinailler. Il y a une coquille dans la formule des float : ce devrait être
(-1)^S, et non-1^S. Cette dernière vaut toujours-1.L’explication selon laquelle l’interprétation du motif de bits brut serait une approximation linéaire par morceaux du logarithme n’est pas exacte non plus. Les lignes entre les points de données du graphe bleu n’existent pas réellement, et il est impossible que seulement la moitié des bits soient à 1. C’est plutôt une version discrète du logarithme ; les points de données qui existent réellement, c’est-à-dire les points où les lignes rouge et bleue se rencontrent, correspondent littéralement à un logarithme mis à l’échelle et translaté. À part cela, c’est un bon article.
[010000, 010111]contient 2, 2,25, 2,5, 2,75, 3, 3,25, 3,5 et 3,75.Mais les mantisses impliquées par le logarithme en base 2 de ces nombres sont respectivement
.0000000,.0010101,.0101001,.0111010,.1001010,.1011001,.1100111et.1110100, et à l’exception de la première, elles ne correspondent pas aux valeurs float001,010, etc. Les float de l’intervalle[2,4)sont espacés linéairement, alors que les logarithmes correspondants ne le sont pas ; on peut donc bien voir les float, comme le dit l’article, comme une approximation linéaire par morceaux du logarithme.Sur le graphe complet, il y aurait
2^32possibilités dans le motif linéaire par morceaux, mais ce n’est pas ce graphe complet que l’article original a dessiné. Comme le texte traite d’entiers 32 bits et d’opérations float IEEE-754 32 bits, je trouve acceptable d’omettre le mot « discret » dans l’explication.C’est un bon article qui explique beaucoup de concepts intéressants, mais le développement algébrique d’une section est étonnamment mauvais.
Après « les étapes exactes pour passer de la première forme à celle-ci sont nombreuses, mais je les ai toutes incluses par souci d’exhaustivité », le développement contient beaucoup d’étapes inutiles, ainsi que plusieurs erreurs de signe qui se compensent. En particulier, le signe moins n’est pas correctement distribué lors du passage de la deuxième à la troisième ligne. À partir de la deuxième ligne, on pourrait partir de
y_n+1 = y_n + (1 - x * y_n^2) / y_n^2 * (y_n^3 / 2)et arriver beaucoup plus rapidement ày_n+1 = y_n (1.5 * y_n - 0.5 * x * y_n * y_n), avec des étapes intermédiaires correctes. Pour quelqu’un qui comprend l’algèbre, ce ne sont, à mon avis, que des étapes évidentes.Le nombre magique du célèbre extrait de code n’est pas la constante optimale. En utilisant une autre constante, on peut probablement réduire l’erreur relative d’environ 0,5 %.
À l’époque, trouver la valeur absolument optimale était peut-être difficile, mais aujourd’hui c’est relativement simple. Je suis moi-même tombé un jour dans ce terrier de lapin et j’ai un notebook Jupyter qui cherche les nombres magiques optimaux pour
(1/x^2)et(1/x).Ce que j’ai trouvé le plus intéressant dans cet article, c’est le lien « How Java's Floating-Point Hurts Everyone Everywhere » : https://people.eecs.berkeley.edu/~wkahan/JAVAhurt.pdf
Son auteur est William Kahan, également connu comme l’« Old Man of Floating-Point » : https://news.ycombinator.com/item?id=29042853 - An Interview with the Old Man of Floating-Point (1998)
JAVAhurtet la mise en page est atroce. On dirait qu’il utilise un package TeX qui étire excessivement, et de façon irrégulière, les espaces entre les mots, ou bien qu’un autre document a été OCRisé avec des espaces supplémentaires insérés.Il y a aussi des espacements supplémentaires bizarres dans les parties en police à chasse fixe. C’était vraiment difficile de rester concentré en le lisant, et même si je sais que ce n’est pas le cas, cela donnait presque l’impression d’un manifeste de geek scientifique.
Cette vidéo que j’avais vue il y a quelque temps était vraiment excellente : https://www.youtube.com/watch?v=p8u_k2LIZyo