Racine carrée inverse rapide

Considérons un nombre x = 0,156 25, pour lequel on souhaite calculer 1/√x ≈ 2,529 82. 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 32-bit :

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 y = 2,614 86 ayant une erreur relative d'environ 3,4 %. Après une itération de la méthode de Newton, le résultat final est y = 2,525 49 avec une erreur de seulement 0,17 %.

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

${\displaystyle {\begin{aligned}x&=\pm 1,b_{1}b_{2}b_{3}\ldots \times 2^{e_{x}}\\&=\pm 2^{e_{x}}(1+m_{x})\end{aligned}}}$

Où l'exposant e_x est un entier, m_x ∈ [0, 1), et 1,b₁b₂b₃... est la représentation binaire de la « mantisse » (1 + m_x). 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 :

• S_x, le « bit de signe », valant 0 si x > 0 ou 1 si x < 0 (1 bit) ;
• E_x = e_x + B est l'« exposant biaisé » où B = 127 est le « biais d'exposant » (en)[12] (8 bits) ;
• M_x = m_x × L, où L = 2²³ [13] (23 bits).

Ces valeurs sont ensuite condensées de gauche à droite dans un conteneur 32 bits.

Par exemple, en utilisant le nombre x = 0,156 25 = 0,001 01₂. En normalisant x on a :

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

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

• S = 0 ;
• E = −3 + 127 = 124 = 01111100₂ ;
• M = 0,25 × 2²³ = 2097152 = 01000000000000000000000₂.

Ces champs sont condensés comme ceci :

S'il fallait calculer 1/√x sans un ordinateur ou une calculatrice, une table de logarithmes serait utilement accompagnée de l'identité log_b(1/√x) = −½ log_b(x) valide quelle que soit la base b. 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 x un nombre normal positif :

${\displaystyle x=2^{e_{x}}(1+m_{x})}$

On a alors

${\displaystyle \log _{2}(x)=e_{x}+\log _{2}(1+m_{x})}$

Mais puisque m_x ∈ [0, 1), le logarithme de la partie droite peut être arrondi par[14]

${\displaystyle \log _{2}(1+m_{x})\approx m_{x}+\sigma }$

où σ est un paramètre arbitraire permettant de régler l'arrondi. Par exemple : σ = 0 fournit des résultats exacts aux bords de l'intervalle tandis que σ ≈ 0.0430357 fournit l'approximation optimale.

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

Alors nous avons l'approximation

${\displaystyle \log _{2}(x)\approx e_{x}+m_{x}+\sigma .}$

D'un autre côté, en interprétant la représentation binaire de x en tant qu'entier, on obtient[15] :

${\displaystyle {\begin{aligned}I_{x}&=E_{x}L+M_{x}\\&=L(e_{x}+B+m_{x})\\&\approx L\log _{2}(x)+L(B-\sigma ).\end{aligned}}}$

On remarque alors que I_x est une approximation linéaire mise à l'échelle et décalée de log₂(x), comme présenté sur le graphique ci-contre. En d'autres termes, log₂(x) est approché par

${\displaystyle \log _{2}(x)\approx {\frac {I_{x}}{L}}-(B-\sigma ).}$

Le calcul de y = 1/√x est fondé sur l'identité

${\displaystyle \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 à x et y, l'équation devient :

${\displaystyle 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

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

à partir duquel on déduit σ ≈ 0,045 046 6. Le second terme, ½ I_x, est déterminé en décalant à droite une fois les bits de I_x[16].

Après avoir appliqué ces opérations, l'algorithme considère de nouveau le mot double comme nombre flottant (y = *(float*)&i;) et effectue une multiplication en nombre flottant (y = y*(1.5f - xhalf*y*y);). 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 :

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

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

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

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

Ainsi y = y*(1.5f - xhalf*y*y); est semblable à ${\displaystyle \,y_{n+1}=y_{n}(1.5-{\frac {xy_{n}^{2}}{2}})={\frac {y_{n}(3-xy_{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 (${\displaystyle y_{n+1}}$) comme argument pour l'itération suivante fait converger l'algorithme sur la racine avec une incertitude de plus en plus faible[17]. Une seule itération a été utilisée dans le cadre du moteur de Quake III, une seconde itération ayant été commentée et laissée.