Le calcul pour un flottant non-dénormal

Pour un flottant positif, on a :

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

Si on prend le logarithme en base 2, on a :

${\displaystyle \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 :

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

Ce qui se simplifie en :

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

Vu que la mantisse est comprise entre 0 et 1, on a : ${\displaystyle \log _{2}(1+m_{x})\approx m_{x}+\sigma }$, avec ${\displaystyle \sigma }$ un terme qui sert à corriger le résultat final. En injectant dans l'expression précédente, on obtient la formule suivante.

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

Avec les nombres flottants, la mantisse m et le terme de correction ${\displaystyle \sigma }$ 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 :

${\displaystyle \log _{2}{x}=e}$

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

L'écriture binaire d'un flottant est égale à une constante près, à un multiple de son logarithme en base 2.

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 ${\displaystyle I}$ dans ce qui suit. Elle est composée de l'exposant placé à gauche de la mantisse. Dit autrement, pour une mantisse de ${\displaystyle n}$ bits, on a :

${\displaystyle I=E<<n+M=E\times 2^{n}+m=E\times L+M}$

On voit que l'écriture binaire est égale au logarithme, mais multiplié par un facteur ${\displaystyle L=2^{n}}$. Le tout tenant dans un entier codé en binaire simple, on a donc bien un nombre encodé en virgule fixe. Cependant, l'exposant ${\displaystyle e}$ est codé en représentation par excès sous la forme ${\displaystyle E=e+{\text{biais}}}$, ce qui donne :

${\displaystyle I=(e+{\text{biais}})\times L+M}$

La mantisse${\displaystyle m}$, au sens mathématique, est censée être comprise entre 0 et 1. Une mantisse de ${\displaystyle n}$ bits est censée coder les ${\displaystyle n}$ rangs après la virgule. Mais dans le nombre flottant, la mantisse ${\displaystyle M}$ occupe les ${\displaystyle n}$ rangs de poids faible avant la virgule. La mantisse ${\displaystyle M}$ contenue dans le nombre flottant est en réalité une version décalée par n rangs de la mantisse ${\displaystyle m}$. Ce décalage est équivalent à une multiplication par ${\displaystyle L=2^{n}}$. On a donc : ${\displaystyle M=2^{n}\times m=L\times m}$. En injectant cette équation dans la précédente, on trouve :

${\displaystyle I=(e+{\text{biais}})\times L+L\times m}$

On factorise :

${\displaystyle I=(e+{\text{biais}}+m)\times L}$

Divisons par ${\displaystyle L=2^{n}}$ :

${\displaystyle {\frac {I}{L}}=e+m+{\text{biais}}}$

On peut alors égaliser cette formule avec l'équation vue plus haut : ${\displaystyle \log _{2}(x)\approx e+m+\sigma }$.

${\displaystyle e+m={\frac {I}{L}}-{\text{biais}}=\log _{2}(x)-\sigma }$

Ce qui donne une approximation du logarithme binaire en virgule fixe, avec un facteur de conversion égal à L.

${\displaystyle \log _{2}(x)={\frac {I}{L}}-{\text{biais}}+\sigma }$

Passage par le logarithme

La première méthode se base sur l'identité mathématique suivante :

${\displaystyle \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.

${\displaystyle \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 : ${\displaystyle 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. 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 ${\displaystyle I}$ d'un flottant. On a, pour rappel :

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

En injectant dans l'équation ${\displaystyle \log _{2}{\left({\frac {1}{\sqrt {x}}}\right)}=-{\frac {1}{2}}\log _{2}{x}}$, on a :

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

On peut alors déduire la représentation entière du résultat en appliquant la formule ${\displaystyle I=(\log _{2}(x)+bias-\sigma )\times L}$. On a donc, en notant ${\displaystyle I_{r}}$ la représentation entière du résultat :

${\displaystyle I_{r}=\left[-{\frac {1}{2}}({\frac {I}{L}}-bias-\sigma )\times L\right]+(bias-\sigma )\times L}$

${\displaystyle I_{r}=-{\frac {1}{2}}I+{\frac {1}{2}}(bias-\sigma )\times L+(bias-\sigma )\times L}$

${\displaystyle I_{r}=-{\frac {1}{2}}I+{\frac {3}{2}}(bias-\sigma )\times L}$

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

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

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