John Carmack est un informaticien américain à l'origine de nombreux jeux vidéos, dont Wolfenstein 3D, Doom et Quake. Ces jeux vidéos sont tous basés sur un moteur 3D pour lequel Carmack a développé ou popularisé de nombreuses techniques innovantes, à une époque où les ordinateurs n'étaient pas épaulés par de puissances cartes graphiques dédiées pour la 3D.

Une opération importante pour ce type de rendu est la normalisation d'un vecteur non nul : on divise un vecteur par sa norme pour n'en garder que la direction. Pour cela, à partir d'un vecteur ${`\\vec u`}$ , on calcule le carré de sa norme à l'aide du produit scalaire ${`\\|\\vec u\\|^2 = \\vec u \\cdot \\vec u`}$ , puis on multiplie chaque composante du vecteur par l'inverse de la racine carrée de ${`\\|\\vec u\\|^2`}$ .

Or, si les additions et multiplications sont peu coûteuses, même pour des flottants, ce n'est pas le cas de la division et plus encore du calcul de la racine carrée. Il s'agit donc d'une fonction à optimiser en priorité.

Lorsque qu'en 2004 est rendu publique le code source du jeu Quake III Arena, on découvre le code suivant (commentaires d'origine) pour calculer l'inverse d'une racine carrée :

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

Cette portion de code, que John Carmack a ainsi popularisé mais n'a pas créé, trouve son origine dans les années 80. On peut lire une esquisse de son histoire ici ou en anglais là.

Dans la suite de cet article, nous allons tenter de la décortiquer.

La partie facile

Le problème de base est simple. Étant donné un réel ${`x > 0`}$ , on veut calculer ${`\\dfrac 1 {\\sqrt x}`}$ , autrement dit on veut trouver une solution à l'équation d'inconnue ${`y`}$ :

&#123;`f(y) = 0`&#125;

pour la fonction ${`f : y \\mapsto \\dfrac 1 {y^2} - x`}$ .

Une façon de résoudre cette équation de manière approchée est d'utiliser la méthode de Newton : en considérant une suite ${`(y_n)`}$ telle que

&#123;`\\forall \\, n \\in \\mathbf N,\\ y_{n + 1&#125; = y_n - \\frac {f(y_n)}{f'(y_n)}`}

alors sous de bonnes conditions (vérifiées ici), la suite ${`(y_n)`}$ converge vers une solution de l'équation ${`f(y) = 0`}$ , autrement dit ici vers ${`\\dfrac 1 {\\sqrt x}`}$ . Or avec la fonction ${`f`}$ que l'on vient de définir, on a :

&#123;`y - \\frac {f(y)&#125;{f'(y)} = \\frac {3 y} 2 - \\frac {x y^3} 2 = y \\times \\Bigl(\\frac 3 2 - \\frac x 2 \\times y \\times y\\Bigr)`}

C'est donc à une étape de la méthode de Newton que correspond la ligne

y = y * (threehalfs - (x2 * y * y));   // 1st iteration

Reste à trouver une bonne valeur ${`y_0`}$ pour avoir une bonne approximation rapidement. C'est tout l'objet des trois lignes mystérieuses que nous allons étudier ensuite :

i = *(long *) &y;                       // evil floating point bit level hacking
i = 0x5f3759df - (i >> 1);              // what the fuck?
y = *(float *) &i;

Et on pourra remarquer que la valeur de départ obtenue avec cette méthode est tellement bonne que la deuxième itération de la méthode de Newton a été commentée. Nous y reviendrons.

La partie diabolique

Notons tout d'abord que ce code manipule des flottants de type float, codés sur 32 bits, ce qui était la norme à l'époque. La répartition est la suivante :

1 bit pour le signe ;
8 bits pour l'exposant auquel on ajoute 127 ;
23 bits pour la mantisse.

Notons que les flottants codés sur 64 bits, courants de nos jours, ont en langage C le type double.

Dans ces trois lignes, la première permet de considérer le codage du flottant y comme… un entier de 32 bits. Elle se comprend ainsi : &y est l'adresse où est stockée la représentation de y. En appliquant (long *), on fait croire à l'ordinateur qu'à cette adresse se trouve un entier (un entier long, pour être précis, sur 32 bits). Ayant maintenant une adresse où est stockée (soi disant) un entier, on récupère la valeur de cet entier avec l'indirection *.

En particulier, pour obtenir l'adresse d'une valeur stockée en mémoire, on utilise &. À l'inverse, à partir d'une adresse, pour obtenir la valeur qui y est stockée, on utilise * qui est agit comme l'opération inverse de la précédente.

Pour résumer, avec l'instruction

i = *(long *)&y;

la variable i contient la valeur de l'entier correspondant à la représentation de y.

Un peu plus tard, après avoir modifié i, nous le retransformons en flottant avec l'instruction

y = *(float *)&i;

Intéressons-nous maintenant à ce qui se passe entre les deux. Pour cela, commençons par étudier les liens entre i et y. Si y représente le flottant non nul ${`2^e (1 + m)`}$ avec ${`0 \\leq m < 1`}$ , alors le mot de 32 bits correspondant à y se décompose ainsi :

où ${`M = 2^{23} m`}$ est un entier. L'entier représenté par i sera donc :

&#123;`2^{23&#125; (127 + e) + M = 2^{23}(127 + e + m)`}

Ainsi, à un flottant ${`x_1 = 2^{e_1} (1 + m_1)`}$ est associé un entier ${`I_1 = 2^{23}(127 + e_1 + m_1)`}$ . On a de plus

&#123;`\\log_2(x_1) = e_1 + \\log_2(1 + m_1)`&#125;

Si l'on définit la fonction

&#123;`g : x \\mapsto \\log_2(1 + x) - x,`&#125;

alors ${`\\log_2(x_1) = e_1 + m_1 + \\delta_1`}$ (où ${`\\delta_1 = g(m_1)`}$ mesure l'erreur entre ${`m_1`}$ et ${`\\log_2(m_1)`}$ ) et donc

&#123;`I_1 = 2^{23&#125;\\bigl(127 + \\log_2(x_1) - \\delta_1\\bigr).`}

Définissons maintenant ${`x_2 = \\dfrac 1 {\\sqrt{x_1}}`}$ que l'on écrit sous la forme ${`x_2 = 2^{e_2}(1 + m_2)`}$ , on a

&#123;`I_2 = 2^{23&#125;\\bigl(127 + \\log_2(x_2) - \\delta_2\\bigr)`}

Mais comme ${`\\log_2(x_2) = - \\frac 1 2 \\log_2(x_1)`}$ et que

&#123;`\\log_2(x_1) = 2^{-23&#125; I_1 - 127 + \\delta_1`}

on en déduit que

&#123;`I_2 = 2^{23&#125;\\Bigl(127 - \\frac 1 2 \\bigl( 2^{-23} I_1 - 127 - \\delta_1 \\bigr) - \\delta_2\\Bigr) = 2^{23} \\Bigl(\\frac 3 2 127 + \\frac 1 2 \\delta_1 - \\delta_2\\Bigr) - \\frac 1 2 I_1`}

Dans cette formule, le ${`- \\frac 1 2 I_1`}$ est traduit exactement par - (i >> 1) (un décalage de bits vers la droite correspondant à une division entière par 2). Comme on veut une approximation, on cherche à remplacer ${`\\delta_1`}$ et ${`\\delta_2`}$ par une constante ${`\\delta`}$ convenablement choisie à l'avance, pour obtenir une expression de la forme

&#123;`I_2 = 2^{23&#125; \\Bigl(\\frac 3 2 127 - \\frac 1 2 \\delta\\Bigr) - \\frac 1 2 I_1.`}

L'étude de la fonction ${`g`}$ montre qu'elle varie sur ${`[0, 1]`}$ entre ${`0`}$ et environ ${`0.0861`}$ , ce qui donne les bornes pour ${`\\delta`}$ .

Pour ${`\\delta = 0`}$ , on a comme constante :

>>> 2**23 * 3 * 127 // 2
1598029824
>>> "{:x}".format(1598029824)
'5f400000'

alors que pour ${`\\delta = 0.0861`}$ , on trouve :

>>> 2**23 * (3 * 127 / 2 - 0.0861 / 2)
1593474390.4256
>>> "{:x}".format(1593474390)
'5efa7d56'

Reste à trouver, entre 0x5efa7d56 et 0x5f400000, la constante qui va minimiser l'erreur dans le calcul de ${`\\frac 1 {\\sqrt x}`}$ . Mais le code initial nous donne peut-être une indication sur ce que l'on doit trouver.

À la recherche du nombre magique

Pour déterminer la meilleure constante à utiliser, le nombre magique, nous allons calculer de façon exhaustive les erreurs relatives maximales en fonction du nombre magique utilisé, et du nombre d'itérations de la méthode de Newton.

La fonction de base que nous allons utiliser est la suivante, où n indique le nombre d'itérations de la méthode de Newton.

float calc(long seed, float x, int n)
{
    float y;
    long i;
    i = *(long *)&x;
    i = seed - (i >> 1); // what the fuck?
    y = *(float *)&i;
    while (n > 0) {
        y = y * (1.5F - (x * 0.5F * y * y));
        --n;
    }
    return fabsf(y);
}

L'erreur relative est calculée ainsi :

float err(float x, float y)
{
    float inv_rac_x = 1.0F / sqrt(x);
    return fabsf((inv_rac_x - y) / inv_rac_x);
}

Pour déterminer l'erreur relative maximale, il suffit de tester toutes les valeurs de mantisses possibles… pour les deux parités d'exposant (puisque lors du décalage vers la droite, le bit de poids faible de l'exposant devient le bit de poids fort de la mantisse), le reste de la valeur de l'exposant (hors sa parité, donc) ne jouant aucun rôle.

On utilise donc la fonction suivante, qui teste toutes les valeurs pour les exposants ${`0`}$ et ${`1`}$ .

On note que pour l'exposant ${`0`}$ (correspondant aux réels dans ${`[1, 2\\mathclose[`}$ ), on utilise un pas de ${`2^{-23}`}$ correspondant à la taille de la mantisse. Pour passer à l'exposant ${`1`}$ (correspondant aux réels dans ${`[2, 4\\mathclose[`}$ ), on doit doubler le pas utilisé. Sinon, celui-ci serait trop petit après le changement d'exposant, et la boucle ne se terminerait pas.

float getmax(long seed,int n)
{
    float m = 0.0F;
    float x, t;
    float step;

    step = pow(2.0, -23.0);

    for (x = 1.0F; x < 2.0F; x += step)
    {
        t = err(x, calc(seed, x, n));
        if (t > m)
        {
            m = t;
        }
    }

    step *= 2.0;

    for (x = 2.0F; x < 4.0F; x += step)
    {
        t = err(x, calc(seed, x, n));
        if (t > m)
        {
            m = t;
        }
    }
    return m;
}

On trouve les valeurs suivantes :

En affinant la recherche, on trouve qu'avec une itération de la méthode de Newton, l'erreur minimale est d'environ ${`0.00175`}$ pour une valeur de 0x5f375a85. Avec une approche plus formelle et moins empirique, Chris Lomont trouve dans l'article Fast Inverse Square Root la valeur de 0x5f375a86 mais la différence d'erreur relative est vraiment négligeable à ce niveau.

Comment expliquer la constante 0x5f3759df utilisée par John Carmack ? Il semblerait qu'elle minimise l'erreur… au bout de deux itérations de la méthode de Newton. Rappelez-vous la deuxième itération présente dans le code initiale, mais commentée car jugée inutile. Cette affirmation, très plausible, est cependant à prendre avec précaution car les fortes « discontinuités » des erreurs relatives montrent que l'on est proche des limites de la précision permise par le codage des flottants en 32 bits.

Seule une itération de la méthode de Newton a été conservée, car elle permettait déjà d'avoir de très bons résultats (une erreur au maximum de 2 pour mille), suffisamment bons pour que le temps d'une seconde itération ne soit pas suffisamment avantageux.

Pour finir, si John Carmack n'a pas utilisé la meilleure constante possible, cela n'a pas dû être très pénalisant, puisque l'erreur relative maximale est seulement passée de ${`0.00175132`}$ à ${`0.00175234`}$ .

L'astuce (qui n'était pas) de John Carmack

Où l'on explique une astuce diabolique concernant les flottants…

La partie facile

La partie diabolique

À la recherche du nombre magique