Racine carrée inverse rapide

Les calculs de lumière (ici dans OpenArena, un fork de Modèle:Lnobr rom) utilisent la racine carrée inverse rapide pour calculer les angles d'incidence et de réflexion.

La racine carrée inverse rapide (en anglais Modèle:Langue, parfois abrégé Fast InvSqrt() ou par la constante 0x5f3759df en hexadécimal) est une méthode pour calculer Modèle:Math, aussi notée Modèle:Math, l'inverse de la racine carrée d'un nombre à virgule flottante à simple précision sur Modèle:Nobr. L'algorithme a probablement été développé chez Silicon Graphics au début des années Modèle:Date. Il a entre autres été utilisé dans le code source de Modèle:Lnobr rom, un jeu vidéo de tir à la première personne sorti en Modèle:DateModèle:Sfn. À l'époque, le principal avantage de cet algorithme était d'éviter des opérations en virgule flottante, coûteuses en puissance de calcul, en préférant des opérations sur entiers. Les racines carrées inverses sont utilisées pour calculer les angles d'incidence et la réflexion pour la lumière et l'ombre en imagerie numérique.

L'algorithme prend en entrée des flottants de Modèle:Nobr non signés et stocke la moitié de cette valeur pour l'utiliser plus tard. Ensuite, il traite le nombre à virgule flottante comme un entier et lui applique un décalage logique à droite d'un bit et le résultat est soustrait à la valeur « magique » 0x5f3759df. Il en résulte une première approximation de la racine carrée inverse du nombre passé en entrée. En considérant de nouveau les bits comme un nombre à virgule flottante et en appliquant au nombre la méthode de Newton, on améliore cette approximation. Bien que n'assurant pas la précision la plus fine possible, le résultat final est une approximation acceptable de la racine carrée inverse d'un nombre à virgule flottante qui s'exécute quatre fois plus vite qu'une division d'un tel nombre.

Motivation

Les normales à une surface sont largement utilisées dans les calculs d'éclairage et d'ombrage, ce qui nécessite le calcul des normes des vecteurs. Ici, on montre un champ de vecteurs normaux à la surface.

Les racines carrées inverses d'un nombre à virgule flottante sont utilisées pour calculer un vecteur normalisé. En synthèse d'image 3D, ces vecteurs normalisés sont utilisés pour déterminer l'éclairage et l'ombrage. Des millions de ces calculs sont ainsi nécessaires chaque seconde. Avant l'apparition de matériel dédié au TCL, ces calculs pouvaient être lents. Ce fut particulièrement le cas lorsque cette technique a été développée au début des années Modèle:Date où les opérations sur les nombres à virgule flottante étaient plus lentes que les calculs sur entiersModèle:Sfn.

Afin de normaliser un vecteur, on détermine la longueur de celui-ci en calculant sa norme euclidienne : la racine carrée de la somme du carré de ses composantes. Après avoir divisé chaque composante par cette longueur, on obtient alors un nouveau vecteur unitaire pointant dans la même direction.

$\Vert 𝒗\Vert =\sqrt{v_1^2+v_2^2+v_3^2}$ est la norme euclidienne du vecteur, de la même manière que l'on calcule une distance dans un espace euclidien.

$\widehat{𝒗}=𝒗/\Vert 𝒗\Vert$ est le vecteur (unitaire) normalisé. Avec $\Vert 𝒗{\Vert }^2$ représentant $v_1^2+v_2^2+v_3^2$,

$\widehat{𝒗}=𝒗/\sqrt{\Vert 𝒗{\Vert }^2}$, liant le vecteur unitaire à la racine carrée inverse des composantes.

Aperçu du code

Le code source suivant est celui de la fonction Q_rsqrt qui implémente la racine carrée inverse rapide dans Modèle:Lnobr rom. On a retiré du code les directives du Modèle:Lnobr (lignes commençant par #), mais on a laissé les commentaires originaux (introduits par //).

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;
}

Afin de déterminer la racine carrée inverse, un programme calculerait une approximation de Modèle:Math puis appliquerait ensuite une méthode numérique afin de peaufiner le résultat jusqu'à atteindre une erreur d'approximation acceptable. Des méthodes d'extraction de racine carrée du début des années Modèle:Date ont permis d'avoir une première approximation depuis une table de correspondance. Cette nouvelle fonction s'est montrée plus efficace que les tables de correspondance et environ quatre fois plus rapide qu'une division flottante classiqueModèle:Sfn. L'algorithme a été conçu selon la norme Modèle:Lnobr pour les nombres à virgule flottante Modèle:Nobr, mais des recherches de Chris Lomont et ensuite Charles McEniry ont montré qu'il pouvait être implémenté en utilisant d'autres spécifications de nombres à virgule flottante.

Le gain de vitesse apporté par le kludge qu'est la racine carrée inverse rapide vient du traitement du mot double contenant le nombre à virgule flottante considéré comme entier qui est ensuite soustrait à une constante spécifique : 0x5f3759df. L'utilité de cette constante n'étant pas claire à première vue, on la considère alors comme un nombre magiqueModèle:SfnModèle:,Modèle:SfnModèle:,Modèle:SfnModèle:,Modèle:Sfn. Après cette soustraction d'entiers et ce décalage à droite, on obtient un mot double qui, lorsqu'il est considéré comme un nombre à virgule flottante, devient une approximation grossière de la racine carrée inverse du nombre entré. Ensuite, une itération de la méthode de Newton est réalisée afin de gagner en précision et le résultat est retourné. L'algorithme génère des résultats raisonnablement précis en utilisant une seule approximation par la méthode de Newton ; toutefois, il reste plus lent que d'utiliser rsqrtss, une instruction du jeu SSE, sortie elle aussi en Modèle:Date sur les processeurs x86.

Exemple pratique

Considérons un nombre Modèle:Math, pour lequel on souhaite calculer Modèle:Math. Voici les premières étapes de l'algorithme :

0011_1110_0010_0000_0000_0000_0000_0000  Trame binaire de x et i
0001_1111_0001_0000_0000_0000_0000_0000  Décalage à droite d'une position : (i >> 1)
0101_1111_0011_0111_0101_1001_1101_1111  Le nombre magique 0x5f3759df
0100_0000_0010_0111_0101_1001_1101_1111  le résultat de 0x5f3759df - (i >> 1)

En utilisant la représentation IEEE 754 sur Modèle:Nobr :

0_01111100_01000000000000000000000  1.25 * 2^-3
0_00111110_00100000000000000000000  1.125 * 2^-65
0_10111110_01101110101100111011111  1.432430... * 2^+63
0_10000000_01001110101100111011111  1.307430... * 2^+1

En réinterprétant la dernière trame binaire en tant que nombre à virgule flottante on obtient l'approximation Modèle:Math ayant une erreur relative d'environ 3,4 %. Après une itération de la méthode de Newton, le résultat final est Modèle:Math avec une erreur de seulement 0,17 %.

Fonctionnement de l'algorithme

L'algorithme calcule Modèle:Math en effectuant les étapes suivantes :

1. Transforme l'argument Modèle:Math en entier afin d'appliquer une approximation de Modèle:Math ;
2. Utilise cet entier pour calculer une approximation de Modèle:Math ;
3. Transforme celui-ci afin de revenir à un nombre flottant afin d'effectuer une approximation de l'exponentielle base-2 ;
4. Affine l'approximation en utilisant une itération de la méthode de Newton.

Représentation en nombre flottant

Puisque cet algorithme s'appuie fortement sur la représentation bit à bit des nombres à virgule flottante simple précision, un aperçu rapide de ce système est fourni ici. Afin d'encoder un nombre réel non nul Modèle:Math en tant que flottant de simple précision, on commence par écrire Modèle:Math comme un nombre binaire en notation scientifique :

$\begin{array}{cc}x& =\pm 1,b_1{b}_2{b}_3\dots \times 2^{e_x}\\ & =\pm 2^{e_x}(1+m_x)\end{array}$

Où l'exposant Modèle:Math est un entier, Modèle:Math, et Modèle:Math est la représentation binaire de la « mantisse » Modèle:Math. Notons qu'il n'est pas nécessaire d'enregistrer le bit avant la virgule dans la mantisse car il vaut toujours 1. Avec ce formalisme, on calcule trois entiers :

• Modèle:Math, le « bit de signe », valant 0 si Modèle:Math ou 1 si Modèle:Math (Modèle:Nobr) ;
• Modèle:Math est l'« exposant biaisé » où Modèle:Math est le Modèle:Lien (Modèle:Nobr) ;
• Modèle:Math, où Modèle:Math (Modèle:Nobr).

Ces valeurs sont ensuite condensées de gauche à droite dans un conteneur Modèle:Nobr.

Par exemple, en utilisant le nombre Modèle:Math. En normalisant Modèle:Math on a :

$x=+2^{-3}(1+0,25)$

Donc, les trois valeurs entières non signées sont :

• Modèle:Math ;
• Modèle:Math ;
• Modèle:Math.

Ces champs sont condensés comme ceci :

Approcher un logarithme en passant à l'entier

S'il fallait calculer Modèle:Math sans un ordinateur ou une calculatrice, une table de logarithmes serait utilement accompagnée de l'identité Modèle:Math valide quelle que soit la base Modèle:Math. La racine carrée inverse rapide repose sur cette dernière ainsi que sur le fait que l'on puisse effectuer un logarithme approximatif d'un nombre en passant d'un float32 à un entier.

Explications :

Soit Modèle:Math un Modèle:Lien positif :

$x=2^{e_x}(1+m_x)$

On a alors

${\log }_2(x)=e_x+{\log }_2(1+m_x)$

Mais puisque Modèle:Math, le logarithme de la partie droite peut être arrondi parModèle:Sfn :

${\log }_2(1+m_x)\approx m_x+\sigma$

où Modèle:Math est un paramètre arbitraire permettant de régler l'arrondi. Par exemple : Modèle:Math fournit des résultats exacts aux bords de l'intervalle tandis que Modèle:Math fournit l'approximation optimale.

L'entier converti en nombre flottant (en bleu), comparé à un logarithme décalé et mis à l'échelle (en gris).

Alors nous avons l'approximation

${\log }_2(x)\approx e_x+m_x+\sigma .$

D'un autre côté, en interprétant la représentation binaire de Modèle:Math en tant qu'entier, on obtient :

$\begin{array}{cc}{I}_x& =E_{x}L+M_x\\ & =L(e_x+B+m_x)\\ & \approx L{\log }_2(x)+L(B-\sigma ).\end{array}$

On remarque alors que Modèle:Math est une approximation linéaire mise à l'échelle et décalée de Modèle:Math, comme présenté sur le graphique ci-contre. En d'autres termes, Modèle:Math est approché par

${\log }_2(x)\approx \frac{I_x}{L}-(B-\sigma ).$

Première approximation du résultat

Le calcul de Modèle:Math est fondé sur l'identité

${\log }_2(y)=-\frac{1}{2}{\log }_2(x)$

En utilisant l'approximation du logarithme telle que précédemment définie et appliquée à Modèle:Math et Modèle:Math, l'équation devient :

$I_y\approx \frac{3}{2}L(B-\sigma )-\frac{1}{2}{I}_x$

Qui s'écrit en code C :

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

Le premier terme étant le nombre magique

$\frac{3}{2}L(B-\sigma )=\text{0x5f3759df}$

à partir duquel on déduit Modèle:Math. Le second terme, Modèle:Math, est déterminé en décalant à droite une fois les bits de Modèle:Math.

Méthode de Newton

Après avoir appliqué ces opérations, l'algorithme considère de nouveau le mot double comme nombre flottant (Modèle:Nobr) et effectue une multiplication en nombre flottant (Modèle:Nobr). Celle-ci étant une itération de la méthode de Newton permettant de trouver des solutions à une équation donnée. Pour ce même exemple :

$y=\frac{1}{\sqrt{x}}$ est la racine carrée inverse, ou encore, en fonction de y :

$f(y)=\frac{1}{y^2}-x=0$.

Avec $y_{n+1}=y_n-\frac{f(y_n)}{f^{\prime }(y_n)}$ représentant l'expression générale de la méthode de Newton avec $y_n$ comme première approximation,

$y_{n+1}=\frac{y_n(3-x{y}_n^2)}{2}$ est l'expression particulière où $f(y)=\frac{1}{y^2}-x$ et $f^{\prime }(y)=\frac{-2}{y^3}$.

Ainsi y = y*(1.5f - xhalf*y*y); est semblable à $y_{n+1}=y_n(1.5-\frac{x{y}_n^2}{2})=\frac{y_n(3-x{y}_n^2)}{2}$

La première approximation est générée en utilisant les opérations en tant qu'entiers puis fournie aux deux dernières lignes de code de la fonction. Des itérations répétées de cet algorithme en utilisant la sortie de la fonction ($y_{n+1}$) comme argument pour l'itération suivante fait converger l'algorithme sur la racine avec une incertitude de plus en plus faibleModèle:Sfn. Une seule itération a été utilisée dans le cadre du [[Id Tech 3|moteur de Modèle:Nobr romains]], une seconde itération ayant été commentée et laissée.

Histoire et enquête

Le code source de Modèle:Lnobr rom a été diffusé après la QuakeCon Modèle:Date, mais des copies de la racine carrée inverse rapide sont apparues sur Usenet et d'autres forums dès Modèle:Date/Modèle:DateModèle:Sfn. Les spéculations à l'époque pointent comme auteur probable de la fonction John Carmack, cofondateur d'id Software, l'entreprise qui a développé les jeux Quake. Mais Carmack dément la chose et suggère que la fonction a été écrite par Terje Mathisen, un programmeur assembleur talentueux qui a aidé les développeurs d'id Software pour optimiser Quake. Mathisen a effectivement écrit une fonction similaire à la fin des années Modèle:Date, mais les auteurs originaux remontent à plus loin dans l'histoire de l'infographie 3D avec l'implémentation faite par Gary Tarolli pour un Modèle:Lien qui serait l'une des premières utilisations connues. Rys Sommefeldt, auteur de l'enquête, finit par conclure que l'algorithme original est l'œuvre de Greg Walsh de Modèle:Lien en collaboration avec Cleve Moler, fondateur de MathWorksModèle:Sfn.

On ne sait pas comment la valeur exacte du nombre magique a été déterminée. Chris Lomont a développé une fonction pour minimiser l'erreur d'approximation en choisissant le nombre magique R dans un intervalle. Il calcule d'abord la constante optimale pour l'étape d'approximation linéaire, il obtient 0x5f37642f qui est proche de 0x5f3759df, mais avec cette nouvelle constante, il obtient une précision moindre après une itération de la méthode de NewtonModèle:Sfn. Il cherche alors une constante optimale même après une ou deux itérations de la méthode de Newton et obtient la valeur 0x5f375a86 qui se révèle plus précise que l'originale, même après chaque étape d'itérationModèle:Sfn. Il conclut alors en se demandant si la valeur originale a été choisie par dérivation ou par essai-erreurModèle:Sfn. Dans la lancée, Lomont indique aussi que la valeur magique pour flottant double précision Modèle:Nobr Modèle:Lnobr est 0x5fe6ec85e7de30da, mais il a été démontré que la valeur exacte était 0x5fe6eb50c7b537a9Modèle:Sfn. Charles McEniry a effectué une optimisation similaire mais plus sophistiquée sur les valeurs probables pour R. Il cherche d'abord par une méthode par force brute et obtient la valeur déterminée par LomontModèle:Sfn. Il a ensuite essayé de rechercher cette valeur par une méthode de dichotomie et obtient alors la valeur utilisée initialement dans la fonction, ce qui conduit McEniry à penser que cette constante a probablement été obtenue par cette méthodeModèle:Sfn.

Voir aussi

Bibliographie

• Modèle:Lien web.
• Modèle:Lien web.
• Modèle:Lien web.
• Modèle:Lien archive.
• Modèle:Lien archive.

Articles connexes

• Extraction de racine carrée

Liens externes

• Modèle:Lien web.