Dans l’histoire de l’informatique, certains morceaux de code sont plus célèbres que d’autres. Ainsi, impossible pour quiconque ayant trainé un peu ses grolles chez les geeks de ne pas être tombé sur le célébrissime hello world.

#include <stdio.h>

int main(void)
{
    printf("hello, worldn");
    return 0;
}

Cette fonction est un classique lorsque l’on souhaite donner un bref aperçu de la syntaxe d’un langage de programmation. Elle est utilisée pour la première fois dans le livre The C Programming Language de Brian W. Kernighan et Dennis Ritchie. Toute aussi célèbre, de nombreux développeurs ont probablement déjà eu, lors d’un entretien, à réécrire la fameuse suite de Fibonacci (ici en OCamL) :

let rec fib n =
  if n < 2 then 1
  else fib (n - 2) + fib (n - 1)

Si la fonction dont nous allons parler aujourd’hui n’est pas aussi célèbre, elle reste une superstar pour quiconque s’est déjà penché sur le développement de jeux vidéo. Il s’agit de la racine carrée inverse introduite dans le moteur du jeu Quake 3 : arena.

Pourquoi est-elle aussi célèbre ?

Calculer l’inverse de la racine carrée d’un nombre est une opération extrêmement utile en modélisation 3D. C’est un calcul auquel on a notamment recours lorsque l’on cherche à calculer des ombres dynamiques. Pour ceux qui ont fait un peu de maths, c’est utilisé pour normaliser un vecteur, qui est l’objet mathématique par lequel tout jeu vidéo représente un objet 3D. Mais calculer l’inverse d’une racine carré suppose de faire deux opérations qui ne sont pas triviales en informatique : une racine carrée et une division. En fait, la division est la seule opération fondamentale qu’un processeur ne sait pas faire nativement. Pour effectuer une division, ce dernier utilise un algorithme qui consiste en une suite d’additions et de soustractions un peu comme l’algorithme d’Euclide que l’on apprend en primaire. Cet algorithme présente un coût qui croît exponentiellement à mesure que l’opération est complexe. Pour l’opération racine carrée, c’est un peu le même problème sauf que l’opération repose sur une suite de multiplications et de divisions qui vont permettre, à chaque pas, de proposer une approximation plus fine de la racine carrée.

Calculer l’inverse d’une racine carrée nécessite donc intuitivement d’utiliser deux algorithmes itératifs qui reposent sur des opération mathématiques couteuses en ressources. Sauf que lorsque l’on joue à un jeu vidéo, il est rare d’avoir très envie d’être coupé dans son action pour laisser le temps au processeur de calculer une racine carrée inverse. L’opération doit donc être extrêmement rapide et c’est en ça que l’algorithme utilisé dans le moteur de Quake 3 est génial. Le voici :

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

Comme vous pouvez le constater, cette fonction ne boucle pas et n’exploite pas la division. Quelques explications s’imposent sur son fonctionnement. À ce stade, je vous prie de m’excuser si les considérations mathématiques paraissent un peu absconses.

L’algorithme exploite la méthode de Newton-Raphson qui est une méthode rapide pour trouver le zéro d’une fonction. Soit la valeur y pour laquelle f(y) = 0. La fonction racine carrée inverse s’exprime comme y=\frac{1}{\sqrt{x}}, ce qui revient à chercher le zéro de la fonction f(y)=\frac{1}{y^2}-x

Il procède par itérations successives pour trouver le point le plus proche de ce zéro en utilisant la formule récurrente:

y_{n+1} = y_n-\frac{f(y_n)}{f'(y_n)}^[1]

Toute la difficulté réside dans le choix de la première approximation y_0. Cette première approximation, elle est donnée par les trois lignes suivantes :

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

Je ne rentrerais pas dans le détail de leur fonctionnement car elles sont aussi ingénieuses que sibyllines pour qui n’a jamais pris de cours d’architecture des ordinateurs et de maths avancées^[2]. Elles utilisent des effets de bord générés par la conversion entre flottant et entier pour simuler un logarithme et une exponentielle en base 2 et un décalage de bit pour la division. Cette astuce permet de donner une approche grossière y_0=\frac{1}{\sqrt{x}}.

Fort de cette première approximation, il ne reste qu’à calculer \frac{f(y_n)}{f'(y_n)} = \frac{\frac{1}{{y_n}^2}-x}{\frac{-2}{{y_n}^3}} = 1.5-\frac{xy{_n}{^2}}{2} et c’est ce que fait la ligne

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

L’explication est un peu velue mais prouve une chose : il a fallu à son auteur de solides connaissances en architecture des ordinateurs et en algèbre pour pouvoir l’écrire.

Le chiffre magique

Une chose qui doit vous interpeler si vous savez développer est l’opération

0x5f3759df - ( i >> 1 );

qui semble sortie de nulle part. Cette ligne permet de faire une division flottante avec… des nombre entiers ! Si le reste de la fonction n’est qu’une application bête et disciplinée de la méthode de Newton-Raphson, c’est dans cette seule ligne que réside tout le génie de l’algorithme. Ce chiffre magique permet de donner une approximation de la racine carrée inversée avec une marge d’erreur initiale de 3, 4%.

Beaucoup de spéculations ont été avancées sur la manière dont l’auteur aurait choisi cette constante et des travaux de recherche ont été menés ensuite pour trouver une valeur qui soit plus optimale. Le mathématicien Chris Lomont analyse, par exemple les résultats en utilisant la constante

0x5f37642f

. Constante qui se révèle, au final, une meilleure approximation initiale mais moins précise après la première itération de Newton-Raphson. Pour faire court : le choix de la constante est très ingénieuse (source, PDF en anglais). La question de savoir si son auteur l’a trouvée par le calcul ou une série d’essais-erreurs reste ouverte.

Qui est l’auteur ?

Beaucoup de personnes ont été soupçonnées d’être l’auteur de cette fonction. John Carmack, le créateur d’Id Software a été le premier à qui ce travail a été attribué. Michael Abrash, un autre talentueux développeur œuvrant sur la série Quake dans les années 90, a également été suspecté du crime. John Carmack répond cependant en 2004 que le code n’est pas de lui et qu’il ne pense pas qu’il soit d’Abrash non plus. Il évoque plutôt le nom de Terje Mathisen. Dans les années 90, Mathisen était un magicien de l’assembleur sur les processeurs x86 et a participé au développement du moteur de Quake.

Dans sa réponse, Terje explique qu’il a effectivement écrit un algorithme similaire à la même période pour une entreprise suédoise à des fins de simulation de comportement de fluides. Il précise cependant que son algorithme était bien plus complexe car l’entreprise avait besoin d’une plus grande précision. Il dit aussi ne pas reconnaitre le code-style général du moteur de Quake 3 : arena et évoque l’idée que l’algorithme puisse être plus ancien et provenir du MIT. En février 1972, en effet, le MIT publie le HACKMEM (en), un rapport technique comportant un ensemble d’algorithmes et de hacks pour le calcul mathématique et la théorie des nombres.

Une dernière piste mène alors à Gary Tarolli, fondateur de 3dfx Interactive, l’un des pionniers des cartes graphiques 3D et créateur de la célèbre Voodoo. Gary Tarolli travaille aujourd’hui pour Nvidia, qui rachète 3dfx Interactive en 2000 pour 120 millions de dollars. Tarolli admet reconnaitre la fonction et avoir travaillé avec des années auparavant mais nie, lui aussi, en être l’auteur.

Alors qui a bien pu écrire ces lignes ? Ben… en fait personne ne sait. La réponse a depuis été oubliée par l’histoire du développement de jeu vidéo.

Conclusion

Si la méthode d’approximation de la racine carrée inverse était brillante lorsqu’elle a été écrite, elle est aujourd’hui obsolète. Créée à une époque où le calcul sur nombre flottant était très coûteux en ressources et où les architectures 32 bits dominaient le marché de l’informatique, elle sera surpassée en 1999 par l’instruction

rsqrtss

introduite par Intel dans l’architecture x86. Aujourd’hui, les calculs en nombre flottant n’ont pratiquement plus aucun coût et le calcul de l’inverse d’une racine carré n’est plus un problème épineux pour les moteurs de jeu.

Mais la fonction

Q_rsqrt

restera tout de même dans l’histoire de l’informatique comme une pierre angulaire du développement de jeu 3D des années 90.

Sources:

Article d’enquête original sur beyon3d.com.
Article Wikipédia sur la fonction (en anglais)
Article Wikipédia sur la méthode de Newton-Raphson.