Les opérations bit à bit/Les opérations arithmétiques sur des flottants

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).

La représentation binaire d'un flottant et les astuces associées

Pour rappel, le codage d'un nombre flottant demande d'écrire celui-ci sous la forme suivante :

$\pm (1+m)\times 2^{e}$, avec $m\in [0,1[$ la mantisse et e l'exposant.

Il ne s'agit ni plus ni moins que l'écriture scientifique que vous avez certainement vu dans votre scolarité, mais adaptée en base 2. En gros, on écrit le nombre sous la forme d'un produit entre une puissance de deux et un autre nombre. L'autre nombre doit être compris entre 1 et 2, pour éviter que l'on ait plusieurs représentations identiques, afin que l'exposant et le nombre multiplié soient uniques. Il est inutile d'encoder la partie entière du nombre, qui vaut toujours 1, mais il est nécessaire d'encoder sa partie fractionnaire. Cela fait donc trois informations à encoder : le signe (flottant positif ou négatif), l'exposant, et la partie fractionnaire du nombre multipliée qui est appelée mantisse.

Le tout est encodé en binaire, en concaténant trois informations : un bit de signe, un champ pour l'exposant et un champ pour la mantisse.

IEEE754 Format Général

Dans ce qui suit, la représentation binaire du flottant sera notée I, celle de l'exposant sera notée E, celle de la mantisse sera notée M. Si on néglige le bit de signe, un nombre flottant se calcule donc comme suit, à partir des deux représentations binaires de l'exposant et de la mantisse :

$I=E<<n+M$

$I=E\times L+M$, avec $L=2^{n}$

La mantisse est codée normalement, sous la forme d'un nombre en virgule fixe/d'un entier. L'exposant est quant à lui un nombre pouvant être nul, positif ou négatif. Mais il n'est pas codé en complément à deux, mais en représentation par excès. La représentation par excès consiste à ajouter un biais aux nombres à encoder afin de les encoder par un entier positif. Pour encoder tous les nombres compris entre -X et Y en représentation par excès, il suffit de prendre la valeur du nombre à encoder, et de lui ajouter un biais égal à X. Ainsi, la valeur -X sera encodée par zéro, et toutes les autres valeurs le seront par un entier positif. Par exemple, prenons des nombres compris entre -127 et 128. On va devoir ajouter un biais égal à 127, ce qui donne :

Valeur avant encodage	Valeur après encodage
-127	0
-126	1
-125	2
…	…
0	127
…	…
127	254
128	255

Pour résumer, partons d'un flottant écrit comme suit : $e=\pm (1+m)\times 2^{e}$. L'exposant est en représentation par excès, ce qui veut dire que l'on a :

$E=e+B$, avec e l'exposant et B le biais de la représentation par excès.

La représentation binaire de la mantisse M s'obtient en prenant la mantisse m (qui est un nombre à virgule compris entre 0 et 1) et en décalant le tout de n rangs vers la gauche. En clair :

$M=m\times L$

On peut alors combiner le tout avec l'équation $I=E\times L+M$, ce qui donne :

$I=(e+B)\times L+m\times L$

$I=(e+B+m)\times L$

Et à partir de cette représentation binaire, certaines opérations deviennent assez simples. Pour cela, il faut cependant traiter le flottant comme s'il s'agissait d'un entier, afin d'utiliser des opérations logiques dessus. Mais en C et dans les autres langages de ce type, les opérations bit à bit ne fonctionnent pas sur les flottants. Il est cependant possible de contourner le problème en faisant quelques manipulations. Dans le langage C, ces manipulations se résument à des conversions et déréférencements entre pointeurs. Le but est alors de copier la représentation binaire d'un flottant dans un nombre entier, sur lequel on pourra effectuer des opérations logiques et/ou mathématiques dessus. Le code pour cela est celui-ci :

    long i  = * ( long * ) &y ; // Bidouille qui permet de copier l'écriture binaire d'un flottant y dans un entier i.

Le code à comprendre est assez simple : on prend un pointeur vers le flottant voulu, puis on convertit ce pointeur sur un flottant en pointeur sur un entier (de type long), et enfin on déréférence ce pointeur.

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. 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 !

(((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.

Et cela peut se comprendre quand on regarde l'équation qui donne la représentation binaire d'un flottant positif :

$I=E\times L+M$, avec $L=2^{n}$

Si on compare deux flottants, le plus grand des deux est celui qui a l'exposant le plus grand, ou la mantisse la plus grande si les deux exposants sont égaux. L'exposant étant dans les bits de poids fort, on comprend que l'exposant passera en priorité dans la comparaison, et que les mantisses seront comparées si les deux exposants sont égaux.

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.

Prenons un nombre x qui vaut :

$x=2^{e}\times (1+m)$, avec m la mantisse et e l'exposant.

Prenons le logarithme en base 2 :

$\log _{2}{x}=\log _{2}{(2^{e}\times (1+m))}$

Vu que le logarithme d'un produit est égal à la somme des logarithmes, on a :

$\log _{2}{x}=\log _{2}{(2^{e})}+\log _{2}{(1+m)}$

Ce qui se simplifie en :

$\log _{2}{x}=e+\log _{2}{(1+m)}$

Le calcul du logarithme d'une somme est complexe et ne peut pas se faire simplement. Heureusement, vu que la mantisse est comprise entre 0 et 1, on a l'approximation suivante qui tient bien la route :

$\log _{2}(1+m)\approx m+\sigma$, avec $\sigma$ un facteur correctif variable, souvent pris égal à 0,0430357.

On a alors :

$\log _{2}{x}\approx e+m+\sigma$

L'équation précédente nous donne une idée de comment calculer le logarithme binaire d'un flottant. Rappelons que les trois termes de l'équation sont des nombres fractionnaires, ce qui fait que les calculs peuvent se faire aussi bien en virgule flottante qu'en virgule fixe, mais pas avec des nombres entiers. Pour calculer le logarithme entier, il faut faire des approximations.

Le logarithme entier d'un nombre flottant

Si l'exposant est codé par un entier dans le nombre flottant, le terme m est un nombre fractionnaire compris entre 0 et 1, de même que $\sigma$. Si on veut un résultat qui soit un nombre entier, on doit se contenter du terme e, ce qui donne :

$\lfloor \log _{2}{x}\rfloor \approx e$

En clair, le logarithme binaire d'un nombre n'est autre que l'exposant de la représentation flottante ! Et pour extraire cet exposant, il suffit de faire quelques décalages sur l'écriture binaire du flottant. Il faut cependant se souvenir que l'exposant est encodé en représentation par excès, ce qui fait qu'il doit lui soustraire le biais adéquat pour obtenir sa valeur normale. La formule mathématique suivante permet de calculer l'exposant à partir de la représentation binaire d'un flottant, notée I, et du biais B de la représentation par excès :

$\lfloor \log _{2}{x}\rfloor \approx E-B$

$\lfloor \log _{2}{x}\rfloor \approx (I>>n)-B$

Le tout peut s'écrire comme ceci :

$\lfloor \log _{2}{x}\rfloor \approx {\frac {I}{L}}-B$, avec $L=2^{n}$.

Précisons qu'ici, la division par L est une division entière.

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

Le calcul du logarithme binaire à partir de l'interprétation de la représentation binaire d'un flottant

Log by aliasing to int

Dans cette section, nous allons voir que le logarithme d'un flottant peut se calculer à partir de sa représentation binaire. Pour simplifier le tout, négligeons totalement le bit de signe pour le moment, en supposant que le nombre traité est positif. Rappelons que l'écriture binaire d'un flottant est mathématiquement définie par la formule suivante, pour un nombre flottant positif :

$I=(e+B+m)\times L$

On peut donc déterminer la somme e + m à partir de la représentation binaire d'un flottant, ce qui est presque égal au logarithme.

${\frac {I}{L}}=e+B+m$

${\frac {I}{L}}-B=e+m$

Mais il reste à faire apparaitre le terme $\sigma$. Pour cela, ajoutons les termes sigma comme ceci :

${\frac {I}{L}}-B+\sigma =e+m+\sigma$

On applique alors la formule $\log _{2}{(x)}\approx e+m+\sigma$ :

${\frac {I}{L}}-B+\sigma \approx \log _{2}{(x)}$

En réarrangeant les termes, on trouve alors :

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

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 : $1 \over {\sqrt {x}}$. L'algorithme utilisé se base sur l'identité mathématique suivante :

$\log {\left({\frac {1}{\sqrt {x}}}\right)}=-{\frac {1}{2}}\log {x}$

Pour simplifier les calculs pour un ordinateur, il est naturel de prendre le logarithme en base 2 de l'opérande.

$\log _{2}{\left({\frac {1}{\sqrt {x}}}\right)}=-{\frac {1}{2}}\log _{2}{x}$

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 : $2^{n}=1<<n$. 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. Avec les flottants, on vient cependant de voir comment calculer le log binaire d'un flottant, ce qui permet de faire la première étape. On vient aussi de voir comment traduire un résultat en virgule fixe en flottant, ce qui correspond à la troisième étape.

Passage par la représentation entière

On vient de voir il y a quelques sections que :

$\log _{2}(x)={\frac {I}{L}}-(B-\sigma )$

On peut combiner les deux équations précédentes pour remplacer le terme de gauche et le terme de droite. Il faut alors faire la différence entre la représentation binaire du résultat, notée $I_{y}$ et celle de l'opérande, notée $I_{x}$. Remplaçons le terme de droite en priorité :

$\log _{2}{\left({\frac {1}{\sqrt {x}}}\right)}=-{\frac {1}{2}}\left[{\frac {I_{x}}{L}}-(B-\sigma )\right]$

Le tout se simplifie en :

$\log _{2}{\left({\frac {1}{\sqrt {x}}}\right)}={\frac {1}{2}}(B-\sigma )-{\frac {1}{2}}{\frac {I_{x}}{L}}$

Maintenant, on peut appliquer la formule $\log _{2}(x)={\frac {I}{L}}-(B-\sigma )$ sur le terme de gauche. On a donc :

${\frac {I_{y}}{L}}-(B-\sigma )={\frac {1}{2}}(B-\sigma )-{\frac {1}{2}}{\frac {I_{x}}{L}}$

On réarrange les termes :

${\frac {I_{y}}{L}}={\frac {3}{2}}(B-\sigma )-{\frac {1}{2}}{\frac {I_{x}}{L}}$

On multiplie par L :

$I_{y}={\frac {3}{2}}(B-\sigma )L-{\frac {1}{2}}I_{x}$

Pour un flottant 32 bits, le terme de gauche vaut : ${\frac {3}{2}}(B-\sigma )\times L={\text{0x5f3759df}}$. On a donc :

$I_{r}={\text{0x5f3759df}}-{\frac {1}{2}}I_{x}$

É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

Imprécision de la racine carré inverse.

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, une méthode utilisée pour déterminer des approximations de fonctions à partir d'une approximation initiale. 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 ) ;