La Racine Carrée Inverse Extrêmement Rapide et Insolite de John Carmack (Quake III)

Question

John Carmack a une fonction spéciale dans le code source de Quake III qui calcule la racine carrée inverse d'un nombre à virgule flottante, 4 fois plus rapidement que la version régulière (float)(1.0/sqrt(x)), incluant une constante étrange 0x5f3759df. Voir le code ci-dessous. Quelqu'un peut-il expliquer ligne par ligne ce qui se passe exactement ici et pourquoi cela fonctionne beaucoup plus rapidement que l'implémentation régulière ?

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

Answer 1

Le code se compose de deux grandes parties. La première partie calcule une approximation pour 1/sqrt(y), et la deuxième prend ce nombre et exécute une itération de la méthode de Newton pour obtenir une meilleure approximation.

Calcul d'une approximation pour 1/sqrt(y)

i  = * ( long * ) &y;
i  = 0x5f3759df - ( i >> 1 );
y  = * ( float * ) &i;

La ligne 1 prend la représentation en virgule flottante de y et la traite comme un entier i. La ligne 2 décale i sur un bit et le soustrait d'une constante mystérieuse. La ligne 3 prend le nombre résultant et le convertit de nouveau en float32 standard. Mais pourquoi cela fonctionne-t-il ?

Soit g une fonction qui associe un nombre en virgule flottante à sa représentation en virgule flottante, lue comme un entier. La ligne 1 ci-dessus définit i = g(y).

Une bonne approximation de g existe(*) :

g(y) ≈ Clog_2 y + D pour certaines constantes C et D. L'intuition pour laquelle une telle bonne approximation existe est que la représentation en virgule flottante de y est approximativement linéaire dans l'exposant.

Le but de la ligne 2 est de passer de g(y) à g(1/sqrt(y)), après quoi la ligne 3 peut utiliser g^-1 pour associer ce nombre à 1/sqrt(y). Utilisant l'approximation ci-dessus, nous avons g(1/sqrt(y)) ≈ Clog_2 (1/sqrt(y)) + D = -C/2 log_2 y + D. Nous pouvons utiliser ces formules pour calculer l'association de g(y) à g(1/sqrt(y)), qui est g(1/sqrt(y)) ≈ 3D/2 - 1/2 * g(y). Dans la ligne 2, nous avons 0x5f3759df ≈ 3D/2, et i >> 1 ≈ 1/2*g(y).

La constante 0x5f3759df est légèrement plus petite que la constante qui donne la meilleure approximation possible pour g(1/sqrt(y)). Cela s'explique par le fait que cette étape ne se fait pas de manière isolée. En raison de la direction dans laquelle la méthode de Newton tend à manquer, l'utilisation d'une constante légèrement plus petite tend à donner de meilleurs résultats. La constante optimale exacte à utiliser dans ce contexte dépend de la distribution d'entrée de y, mais 0x5f3759df est l'une des constantes qui donne de bons résultats sur une plage assez large.

Une description plus détaillée de ce processus peut être trouvée sur Wikipedia : https://fr.wikipedia.org/wiki/Fast_inverse_square_root#Algorithm

(*) Plus explicitement, soit y = 2^e*(1+f). En prenant le log des deux côtés, nous obtenons log_2 y = e + log_2(1+f), qui peut être approximé comme log_2 y ≈ e + f + σ pour une petite constante sigma. Séparément, le codage float32 de y exprimé comme un entier est g(y) ≈ 2^23 * (e+127) + f * 2^23. En combinant les deux équations, nous obtenons g(y) ≈ 2^23 * log_2 y + 2^23 * (127 - σ).

Utilisation de la méthode de Newton

y  = y * ( threehalfs - ( x2 * y * y ) );

Considérons la fonction f(y) = 1/y^2 - num. Le zéro positif de f est y = 1/sqrt(num), ce que nous cherchons à calculer.

La méthode de Newton est un algorithme itératif pour prendre une approximation y_n du zéro d'une fonction f, et calculer une meilleure approximation y_n+1, en utilisant l'équation suivante : y_n+1 = y_n - f(y_n)/f'(y_n).

Calculant ce que cela donne pour notre fonction f donne l'équation suivante : y_n+1 = y_n - (-y_n+y_n^3*num)/2 = y_n * (3/2 - num/2 * y_n * y_n). C'est exactement ce que fait la ligne de code ci-dessus.

Vous pouvez en apprendre davantage sur les détails de la méthode de Newton ici : https://fr.wikipedia.org/wiki/M%C3%A9thode_de_Newton