Après avoir vu comment certains calculs entiers peuvent se réaliser via des opérations logiques, il est temps de voir des astuces de calcul sur des nombres flottants. Dans ce qui va suivre, nous allons travailler avec des nombres flottants au format IEEE754, le format le plus courant (pour ne pas dire le seul utilisé à l'heure actuelle).
Le calcul de la valeur absolue d'un flottant.
Nous allons commencer par une opération extrêmement simple : la valeur absolue. On sait que les flottants sont codés avec une représentation spéciale. Pour rappel, un nombre flottant encode est encodé avec un bit de signe, un champ pour l'exposant et un champ pour la mantisse. La valeur décimale du flottant se calcule comme suit :
Comme on le voit, le signe d'un flottant est indiqué par un bit de signe, localisé dans le bit de poids fort. Obtenir sa valeur absolue demande simplement de mettre ce bit de signe à 0. Un simple masque suffit ! Cependant, cela demande de travailler sur l'écriture binaire du flottant, chose qui n'est pas possible sans quelques tricks assez velus. En C et dans les autres langages de ce type, les opérations bit à bit ne fonctionnent pas sur les flottants, sans compter que le masque est obligatoirement un entier. Pour travailler sur l'écriture binaire, le mieux est de passer par des conversions de pointeurs : on prend un pointeur sur le flottant à traiter, on le convertit en pointeur vers un unsigned/int, et on déréférence : ob a alors l'écriture binaire du flottant directement dans un entier (si le flottant et l'entier ont des tailles identiques). Le code est donc le suivant :
(((int *) &x) + 1) &= 0x7fffffff ;
Les comparaisons entre flottants.
Chose peu connue, il est possible de comparer deux nombres flottants en utilisant les comparaisons pour nombres entiers. Du fait de leur encodage, les comparaisons entre flottants et entiers sont différentes sur le papier. Mais dans les faits, quelques propriétés de l'encodage des flottants fait que l'on peut comparer deux flottants comme s'ils s'agissait d'entiers sans problèmes. Ce que cela veut dire, c'est que si l'on prend l'écriture binaire de deux flottants, et que l'on compare ces écritures comme s'il s'agissait d'entiers, alors le résultat de la comparaison sera le bon. Cela vient du fait que deux flottants consécutifs auront un encodage consécutifs, chose qui est garantie par la norme IEEE754. Tout du moins, cela fonctionne dans la plupart des cas, bien que cela demande quand même quelques précisions. Premièrement, l'encodage du bit de signe pose quelques problèmes : les nombres négatifs suivent les positifs du point de vue de l'écriture binaire. La solution est alors d'utiliser une comparaison signée. Pour régler ce problème, il suffit d'inverser le bit de signe avec le masque adéquat.
Le calcul du logarithme d'un flottant.
Nous verrons dans quelques chapitres comment trouver le logarithme d'un entier, pour des raisons liée à la nature de ce calcul, mais la méthode que nous aborderons ne fonctionne évidemment pas sur les nombres flottants. Pour ceux-ci, il existe cependant une méthode assez simple, qui demande de passer par l'écriture binaire du flottant en question. La technique que l'on va voir donne un résultat entier, pas un flottant.
Le calcul pour un flottant non-dénormal.
Pour un flottant positif, on a :
Si on prend le logarithme en base 2, on a :
Vu que le logarithme d'un produit est égal à la somme des logarithmes, on a :
Ce qui se simplifie en :
Vu que la mantisse est comprise entre 0 et 1, on a : formula_5, avec formula_6 un terme qui sert à corriger le résultat final. En injectant dans l'expression précédente, on obtient la formule suivante.
Avec les nombres flottants, la mantisse m et le terme de correction formula_6 sont des nombres fractionnaires, compris entre 0 et 1, ce qui fait qu'on les éliminera dans les calculs qui vont suivre. Mais si on travaillait en virgule fixe, nous devrions les prendre en compte. Nous le ferons d'ailleurs dans la suite de ce cours. En négligeant ces termes, on trouve alors :
Pour faire ce calcul, il suffit de prendre l'écriture binaire du flottant, de faire quelques décalages. Il faut cependant se souvenir que l'exposant est encodé en représentation par excès, ce qui fait qu'il doit lui soustraire un biais pour obtenir sa valeur normale. Le code devrait être le suivant pour un flottant 32 bits (simple précision). La première ligne permet de traiter directement l'écriture binaire du flottant, afin de faire les décalages et autres. Il faut cependant signaler que ce code ne fonctionne pas avec les flottants dénormaux.
float rsqrt(float n)
 long i = * ( long * ) &n ;
 unsigned e = (c » 23) - 127 ;
 return e ;
L'interprétation de l'écriture binaire d'un flottant.
Nous allons terminer cette section par une petite anecdote, qui sera utile dans la suite du cours. Nous allons voir que la représentation binaire d'un flottant a un lien très profond avec son logarithme. Pour faire simple, l'écriture binaire du flottant approxime très bien le logarithme binaire encodé en virgule fixe. Dit autrement, l'écriture binaire d'un flottant est égale à une constante près, à un multiple de son logarithme en base 2. La représentation binaire du flottant sera notée formula_10 dans ce qui suit. Elle est composée de l'exposant placé à gauche de la mantisse. Dit autrement, pour une mantisse de formula_11 bits, on a :
On voit que l'écriture binaire est égale au logarithme, mais multiplié par un facteur formula_13. Le tout tenant dans un entier codé en binaire simple, on a donc bien un nombre encodé en virgule fixe. Cependant, l'exposant formula_14 est codé en représentation par excès sous la forme formula_15, ce qui donne :
La mantisseformula_17, au sens mathématique, est censée être comprise entre 0 et 1. Une mantisse de formula_11 bits est censée coder les formula_11 rangs après la virgule. Mais dans le nombre flottant, la mantisse formula_20 occupe les formula_11 rangs de poids faible avant la virgule. La mantisse formula_20 contenue dans le nombre flottant est en réalité une version décalée par n rangs de la mantisse formula_17. Ce décalage est équivalent à une multiplication par formula_24. On a donc : formula_25. En injectant cette équation dans la précédente, on trouve :
On factorise :
Divisons par formula_13 :
On peut alors égaliser cette formule avec l'équation vue plus haut : formula_7.
Ce qui donne une approximation du logarithme binaire en virgule fixe, avec un facteur de conversion égal à L.
Le calcul de la racine carrée inverse d'un flottant.
Dans cette section nous allons nous intéresser au calcul de la racine carrée inverse, à savoir : formula_33. L'algorithme utilisé pour la calculer s'appelle la méthode de Newton/Raphson, ou encore méthode de Newton. Sans rentrer dans des considérations mathématiques trop compliquée, celle-ci a besoin d'une approximation pas trop mauvaise du résultat, qui sera notée formula_34. À partir de cette approximation, cette méthode applique un calcul de manière répétée, le résultat du calcul s'approchant de plus en plus du résultat réel. Cette méthode marche pour de nombreuses fonctions, tel le logarithme, l'exponentielle, et bien d'autres, et est naturellement utilisée pour la racine carrée inverse. Dans le cas de la racine carrée inverse, ce calcul est le suivant :
Reste à calculer l'approximation de base formula_34. Pour cela, il existe diverses méthodes que nous allons maintenant aborder.
Passage par le logarithme.
La première méthode se base sur l'identité mathématique suivante :
Pour simplifier les calculs pour un ordinateur, il est naturel de prendre le logarithme en base 2 de l'opérande.
La solution consiste donc à calculer le logarithme de l'opérande, multiplier par - 0.5, et prendre le logarithme inverse (l'exponentielle). On sait déjà comment calculer ce logarithme en base 2 pour un entier, la division ne pose pas de problèmes (c'est un simple décalage) et le calcul de l'exponentielle est trivial : formula_39. Mais cette méthode ne fonctionne que pour des entiers : il faut donc convertir le nombre flottant voulu en entier et appliquer dessus les calculs. Ce qui pose de gros problèmes : il faudrait des entiers extrêmement longs, de plusieurs centaines de bits, pour pouvoir faire le calcul sans approximation extrême. La solution qui consiste à faire les calculs sur les flottants est encore pire, vu que le calcul du logarithme et de l'exponentielle se font alors par la méthode de Newton et sont donc aussi rapides/lents que le calcul de la racine carrée inverse. Mais comme nous allons le voir dans la suite, il existe un moyen de ruser !
Passage par la représentation entière.
Pour cela, il faut juste réutiliser ce qu'on a vu dans la section sur le logarithme. On sait qu'on peut donner une approximation du logarithme à partir de la représentation entière formula_10 d'un flottant. On a, pour rappel :
En injectant dans l'équation formula_38, on a :
On peut alors déduire la représentation entière du résultat en appliquant la formule formula_44. On a donc, en notant formula_45 la représentation entière du résultat :
Pour un flottant 32 bits, le terme de droite vaut : :formula_49. On a donc :
Évidemment, il ne s'agit là que d'une approximation. Le code correspondant est le suivant :
float rsqrt(float number)
 long i = * ( long * ) &y ; // Bidouille qui permet de traiter l'écriture binaire d'un flottant.
 i = 0x5f3759df - ( i » 1 ) ; // Utilisation du calcul précédent, en remplaçant la division par deux par un décalage.
 float y = * ( float * ) &i; // Conversion du résultat en flottant, facultatif.
 return y;
Le calcul de la racine carrée inverse rapide d'un flottant.
Les raisonnements précédents sont à 'origine d'une bidouille vraiment intéressante, incorporée dans certains moteurs de jeux vidéo. Celle-ci se base sur le principe vu dans la section précédente. On ne sait pas qui l'a inventé, mais il est devenu célèbre lors de la publication du code source de Quake 3 arena, un vieux fast FPS comme on en fait plus, très sympa à jouer en LAN. Un vrai jeu, quoi ! Bref, dans les jeux vidéo, on a souvent besoin de calculer la racine carrée inverse. Ce calcul est évidemment effectué sur des nombres flottants. Mais à l'époque de la création de ce jeu vidéo, les opérations flottante étaient très lentes. Et elles le sont toujours un petit peu, d'ailleurs. Pour diminuer la lenteur du calcul de la racine carrée inverse, on a inventé la Fast Inverse SQRT, afin de calculer celle-ci en utilisant uniquement des opérations entières ! Cet algorithme prend un flottant 32 bits, et donne une approximation de la racine carrée inverse. La qualité de cette approximation est illustrée à droite. Cet algorithme est très simple, jugez plutôt :
float 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;
La version en C moderne, à laquelle on aurait omis quelques optimisations, donnerait ceci :
float rsqrt(float number)
 float y = number;
 long i = * ( long * ) &y; // evil floating point bit level hacking
 i = 0x5f3759df - ( i » 1 ); // what the fuck? 
 y = * ( float * ) &i;
 y = y * ( 1.5F - ( y * y * number / 2 ) ); // 1st iteration
 //y = y * ( 1.5F - ( y * y * number / 2 ) ); // 2nd iteration, this can be removed
 return y;
Les deux dernières lignes de code, dont celle mise en commentaire, correspondent chacune à une itération de la méthode de Newton. On peut remarquer que l'auteur de cette routine a appliqué quelques optimisations sur ce calcul : il a notamment remplacé la division par 2 par une multiplication par 0.5, les multiplications étant plus rapides que les divisions. Les lignes de code précédentes calculent une approximation assez précise de la racine carrée inverse, qui est ensuite utilisée telle quelle par la méthode de Newton. Le premier "bloc" de code contient de simples déclarations de variables utilisées dans le calcul, la plupart servant uniquement à rendre les calculs plus lisibles. Les lignes suivantes contiennent des conversions et déréférencements entre pointeurs qui permettent d’accéder directement à l'écriture binaire du flottant, afin de traiter celle-ci avec des opérations logiques comme s'il s'agissait d'un entier. Le calcul de l'approximation réutilise ce qu'on a vu dans la section précédente, par la ligne de code suivante :
i = 0x5f3759df - ( i » 1 ) ;