Maîtrisez les nombres à virgule en C

La base de la base

Les types pour représenter un nombre à virgule

1
2

printf("taille d'un float : %u\n", sizeof(float));
printf("taille d'un double : %u\n", sizeof(double));

Intérêt des nombres flottants

Écrire une constante

1
2
3
4

0.0e0 ;     0.0 ;
 0.e0 ;      0. ;
 .0e0 ;      .0 ;
  0e0 ;

• la lettre 'F' minuscule ou majuscule demande un float ;
  
• la lettre 'L' minuscule ou majuscule demande un long double.
  

• ajoutez le préfixe 0x (ou 0X) avant votre constante (mais après le signe) ;
• écrivez la partie significative en base 16 (c'est-à-dire avec les chiffres de 0 à 9 et les lettres de A à F) ;
• l'exposant est obligatoire (à la limite, vous n'avez qu'à mettre 0) ; il doit être écrit en base 10 comme d'habitude (et non en hexadécimal), mais précédé de la lettre p (ou P) au lieu de e. De plus, il se réfère à une puissance de 2 et non de 10.

Des constantes inexactes

Quelques informations pratiques

Les fonctions d'entrée et de sortie formatées : printf & scanf

Les fonctions mathématiques avec <math.h>

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105

#include <math.h>


/** fonctions courantes **/

double pow(double x, double n);   // x à la puissance n

double sqrt(double);   // racine carrée
double cbrt(double);   // [C99] racine cubique

double fabs(double);           // valeur absolue
double fdim(double, double);   // [C99] différence positive
double fmin(double, double);   // [C99] minimum
double fmax(double, double);   // [C99] maximum


/** fonctions trigonométriques (les angles sont en radians) **/

double cos(double);   // cosinus
double sin(double);   // sinus
double tan(double);   // tangente
double acos(double);  // arc-cosinus
double asin(double);  // arc-sinus
double atan(double);  // arc-tangente
double atan2(double y, double x);  /* arc-tangente avec 2 arguments : angle entre l'axe
                                      des abscisses et le point de coordonnées (x ; y) ;
                                      attention, y est le 1er argument ! */
double cosh(double);  // cosinus hyperbolique
double sinh(double);  // sinus hyperbolique
double tanh(double);  // tangente hyperbolique
double acosh(double);  // [C99] arc-cosinus hyperbolique
double asinh(double);  // [C99] arc-sinus hyperbolique
double atanh(double);  // [C99] arc-tangente hyperbolique


/** arrondis, parties entière et fractionnaire **/

double          round(double);   // [C99] arrondi à l'entier, tel que round(.5)==1
long int       lround(double);   // [C99] idem
long long int llround(double);   // [C99] idem

double          rint(double);   // [C99] arrondi à l'entier suivant le mode d'arrondi  [1]
long int       lrint(double);   // [C99] idem
long long int llrint(double);   // [C99] idem

double nearbyint(double)   /* [C99] arrondi à l'entier suivant le mode d'arrondi,
                              sans lever d'exception "inexact"  [1] */

double ceil(double);   // arrondi à l'entier supérieur
double floor(double);  // arrondi à l'entier inférieur

double trunc(double);   // [C99] partie entière (troncature)
double modf(double, double*);   /* retourne la partie fractionnaire et écrit la partie
                                   entière dans le nombre pointé */


/** exponentielles, logarithmes & co. **/

double exp(double);       // fonction exponentielle
double expm1(double x);   // [C99] équivalent à exp(x)-1, mais plus précis
double exp2(double x);    // [C99] 2 à la puissance x

double log(double);       // logarithme naturel (ou "népérien")
double log1p(double x);   // [C99] équivalent à log(1+x), mais plus précis
double log10(double);     // logarithme en base 10
double log2(double);      // [C99] logarithme en base 2

double logb(double x);    /* [C99] logarithme entier en base 2 (équivalent à
                             floor(log2(x)), mais plus rapide)  [2] */
int    ilogb(double);     // [C99] idem

double scalbn(double x, int n);     // [C99] x * 2^n  [2]
double scalbln(double, long int);   // [C99] idem
double ldexp(double x, int n);   // x * 2^n
double frexp(double x, int*);    /* fraction normalisée : écrit un exposant EXP dans l'entier
                                    pointé et renvoie une valeur A, tels que x == A / (2^EXP) ;
                                    A est compris entre 0.5 inclus et 1 exclus */


/** plus de fonctions **/

double fmod(double, double);        // reste de la division
double remainder(double, double);   // [C99] idem (avec une subtile nuance)

double hypot(double, double);   /* [C99] hypoténuse du triangle rectangle dont les deux
                                   autres côtés sont donnés */
double erf(double);    // [C99] fonction d'erreur de Gauss
double erfc(double);   // [C99] fonction d'erreur complémentaire
double tgamma(double);   // [C99] fonction gamma
double lgamma(double);   // [C99] logarithme naturel de la valeur absolue de la fonction gamma

double copysign(double x, double y);   // [C99] valeur de x avec le signe de y
double fma(double x, double y, double z)   // [C99] x*y + z

double nextafter(double x, double y);     /* [C99] prochaine valeur flottante représentable
                                             après x (vers y) */
double nexttoward(double, long double);   // [C99] idem


/** NOTES :
[1] Modes d'arrondis, exceptions : on verra ça plus tard.
[2] En fait, j'ai mis 2 pour simplifier, mais en réalité il s'agit de FLT_RADIX
    (une constante définie dans le header <float.h> et dont vous n'avez pas encore
    les connaissances requises pour en comprendre le sens) ; FLT_RADIX vaut probablement 2.
**/

• une qui travaille avec des double : celle présentée ;
• une qui travaille avec des float : il faut ajouter la lettre 'f' à la fin du nom de la fonction ;
• une qui travaille avec des long double : il faut ajouter la lettre 'l' à la fin du nom de la fonction.

• Project > Build options… > onglet Linker settings ;
• dans le champ Link libraries (à gauche), cliquez sur Add et entrez le nom de la bibliothèque à lier, soit « m », puis validez ;
• autre possibilité pour cette dernière étape : taper simplement « -lm » dans le champ Other linker options (à gauche).
• Voilà, les fonctions de <math.h> seront désormais incluses dans vos programmes !

• la manière de représenter les nombres flottants en mémoire avec plusieurs formats ;
• cinq opérations associées : l'addition (+), la soustraction (-), la multiplication (*), la division (/) et la racine carrée (sqrt()) ;
• des modes d'arrondi ;
• des exceptions.

Principes généraux

• le signe, codé sur un bit : ce bit vaut 0 si le nombre est positif, 1 s'il est négatif.
  Remarquez ici que la représentation du signe des nombres flottants se fait selon le bit de poids fort et non selon la règle du complément à 2. En conséquence, il existe deux représentations pour zéro, un positif et un négatif (c'est néanmoins un défaut mineur, les implémentations gérant cette dualité ; par exemple, +0.0==-0.0 renverra vrai), et les valeurs possibles sont totalement "symétriques" par rapport à zéro (c'est la raison pour laquelle j'utiliserai souvent le symbole ±, plus ou moins, par la suite) ;
   
• l'exposant, c'est-à-dire la puissance à laquelle il faut élever 2 (et non 10 !).
  Le nombre de bits qu'il occupe dépend de la taille du type.
  Ce nombre est un entier et peut être négatif. Pour le représenter, on n'utilise ni la règle du bit de signe, ni celle du complément à 2 (car cela rendrait la comparaison de nombres plus difficile, comme on le verra ensuite) ; on opère un décalage de l'exposant réel en lui ajoutant 2^e-1-1, où e représente le nombre de bits occupés par l'exposant ; en gros, ce nombre représente la valeur où tous les bits sauf celui de poids fort sont à 1 (c'est-à-dire 127 = 0b01111111 si l'exposant est codé sur 8 bits, ou 1023 = 0b01111111111 s'il est codé sur 11 bits) ;
  
  Remarque : Par la suite, j'emploierai les termes d'exposant "décalé" et "non-décalé" pour différencier l'exposant tel qu'il est représenté en mémoire de l'exposant ainsi représenté.
   
• la mantisse, c'est-à-dire la valeur de la partie décimale en notation binaire scientifique (avec un seul chiffre avant la virgule). On ne garde que la partie décimale car, du fait de la notation binaire, la partie entière est forcément égale à 1 (sauf dans les cas particuliers de zéro et de certains nombres très petits dont on reparlera). On parle de bit implicite (un bit "caché" représentant la partie entière, qui ici vaudrait 1).
  On économise ainsi un bit, ce qui permet d'avoir une précision plus grande.
  La mantisse occupe les bits de poids faibles restants.

Les différents types de nombres représentables

• l'infini positif ou négatif (noté par les mathématiciens). C'est par exemple le résultat de la division d'un nombre non nul par zéro (le signe dépendant alors du signe du numérateur et du zéro en dénominateur) ;
• NaN (not a number, pas un nombre). NaN est une valeur spéciale qui est utilisée pour signaler une erreur ( ou par exemples). N'importe quel calcul avec un NaN doit renvoyer NaN (sauf quelques exceptions), et n'importe quelle comparaison avec un NaN doit renvoyer faux (sauf !=) ; un NaN n'est même pas égal à lui-même, c'est pourquoi x==x peut renvoyer faux ;
• enfin, on n'a pas encore réglé le problème de zéro.

• Si l'exposant vaut sa valeur maximale (soit 2^e-1) et que la mantisse est nulle, alors c'est l'infini (positif ou négatif, selon le bit de signe).
  
• Si l'exposant vaut sa valeur maximale (soit 2^e-1) et que la mantisse n'est pas nulle, alors c'est NaN.
  
• Si l'exposant décalé vaut 0 et que la mantisse est nulle, alors c'est zéro (positif ou négatif, selon le bit de signe).
  
• Si l'exposant décalé vaut 0 et que la mantisse n'est pas nulle, alors c'est un nombre dénormalisé : on considère que le bit implicite (la partie entière en notation scientifique) vaut 0.
  
• Dans tous les autres cas, c'est un nombre normalisé : on considère que le bit implicite vaut 1.
  

Tableau récapitulatif des différents types de valeurs d'un nombre à virgule flottante

Type	Représentation mémoire		Valeur
	binaire	hexadécimal
Not a Number	[0/1] — 11111111 — 11111111111111111111111 … [0/1] — 11111111 — 00000000000000000000001	[7/F]F FF FF FF … [7/F]F 80 00 01	NaN … NaN
Infini	[0/1] — 11111111 — 00000000000000000000000	[7/F]F 80 00 00	±
Nombre normalisé	[0/1] — 11111110 — 11111111111111111111111 … [0/1] — 00000001 — 00000000000000000000000	[7/F]F 7F FF FF … [0/8]0 80 00 00	± 3,4028235 ×10+38 … ± 1,1754944 ×10-38
Nombre dénormalisé	[0/1] — 00000000 — 11111111111111111111111 … [0/1] — 00000000 — 00000000000000000000001	[0/8]0 7F FF FF … [0/8]0 00 00 01	± 1,1754942 ×10-38 … ± 1,4012985 ×10-45
Zéro	[0/1] — 00000000 — 00000000000000000000000	[0/8]0 00 00 00	± 0

idem pour un double

Type	Représentation mémoire	Valeur
Not a Number	[7/F]F FF FF FF   FF FF FF FF … [7/F]F F0 00 00   00 00 00 01	NaN … NaN
Infini	[7/F]F F0 00 00   00 00 00 00	±
Nombre normalisé	[7/F]F EF FF FF   FF FF FF FF … [0/8]0 10 00 00   00 00 00 00	± 1,7976931348623157 ×10+308 … ± 2,2250738585072014 ×10-308
Nombre dénormalisé	[0/8]0 0F FF FF   FF FF FF FF … [0/8]0 00 00 00   00 00 00 01	± 2,2250738585072010 ×10-308 … ± 5,0000000000000000 ×10-324
Zéro	[0/8]0 00 00 00   00 00 00 00	± 0

Nombres dénormalisés

Mais pourquoi avoir inventé les nombres dénormalisés ? Il aurait été plus simple de décider qu'à l'exposant décalé nul correspondaient également des nombres normalisés "banals" (dont l'exposant réel serait bien -127 et non -126).

• Historiquement, l'exposant minimal servait pour représenter zéro, quelle que soit la valeur de la mantisse ; les nombres dénormalisés ont été introduits en 1985 avec la première version de la norme IEEE 754.
• D'ailleurs, si l'on avait choisi des nombres normalisés pour "remplir" cet exposant, la représentation de zéro n'aurait plus été cohérente. En effet, on se serait retrouvé avec ±zéro = ± 1,00000000000000000000000 ×2^-127 (le bit implicite valant 1).
  Cela dit, les nombres dénormalisés introduisent eux-même une incohérence, puisqu'il y a deux types de nombres qu'il faut interpréter différemment ; de plus, on aurait été tout à fait capable de gérer cette incohérence pour interpréter correctement zéro.
• Les nombres dénormalisés permettent de représenter des nombres bien plus petits (plus proches de zéro) que si l'on avait choisi des nombres normalisés. En effet, le plus petit nombre dénormalisé est ± 0b 0,00000000000000000000001 ×2^-126, soit ±1 ×2^-149 (-126-23 = -149) ; si l'on avait choisi des nombres normalisés, le plus petit aurait été ± 0b 1,00000000000000000000001 ×2^-127.

Retrouver la valeur du nombre flottant

• la valeur du flottant représenté ;
• et le nombre de bits occupés par la mantisse et l'exposant (respectivement) ;
• la mantisse (plus précisément, sa représentation entière) ;
• l'exposant non-décalé (en termes de puissances de 2) ; le décalage de l'exposant est , on en déduit les valeurs limites de l'exposant :
  • la valeur minimale de l'exposant non-décalé, réservée aux zéros et nombres dénormalisés, est (mais l'exposant "réel" des dénormalisés est ) ;
  • la valeur maximale de l'exposant non-décalé, réservée aux infinis et NaNs, est .

• 1 ou 0 représente la partie entière du flottant (c'est-à-dire le bit implicite).
• La fraction qui suit représente la partie décimale (on a ).

• (normalisé) ;
• (dénormalisé).

Exemple 1 : Nombre normalisé

Exemple 2 : Nombre dénormalisé

Intervalle entre les nombres en fonction de l'exposant

• 2^-126 - 23 = 2^-149 ≈ 1,40 ×10^-45 : écart entre les nombres dénormalisés et les premiers nombres normalisés (et valeur du tout premier nombre dénormalisé) ;
• 2^0 - 23 ≈ 1,19 ×10^-7 : écart entre les nombres normalisés compris entre 1 et 2 ;
• 1 : écart entre les nombres normalisés d'exposant non-décalé 23 ;
• 2^127 - 23 ≈ 2,03 ×10³¹ : écart entre les plus grands nombres normalisés.

Précision et chiffres significatifs

• le format 32 bits garantit 7 chiffres significatifs en base 10 ;
• le format 64 bits garantit 16 chiffres significatifs en base 10.

un nombre normalisé a donc chiffres significatifs en binaire

Les exceptions

• invalid operation (opération invalide) : se produit lorsqu'une opération interdite est effectuée ; le résultat du calcul en question sera NaN ;
• division by zero (division par zéro) : se produit lorsqu'on tente de diviser un nombre (non nul) par zéro ; le résultat sera ±  ;
• overflow : se produit lorsque le résultat d'un calcul est trop grand (en valeur absolue) pour être stocké ; on arrondit à ±  ;
• underflow : se produit lorsque le résultat d'un calcul est trop petit (en valeur absolue) pour être stocké ; on arrondit à zéro ;
• inexact : se produit lorsqu'on effectue un autre arrondi pour pouvoir stocker un nombre flottant (parce que le résultat a plus de chiffres significatifs qu'on peut en stocker).

Les arrondis

• au nombre le plus proche : c'est le comportement par défaut (si le nombre à arrondir tombe pile-poile entre deux, alors on arrondit à celui dont le dernier bit vaut 0) ;
• vers zéro ;
• vers + ;
• vers -.

1
2
3
4
5
6

#include <stdio.h>
int main(void) {
   printf("%g\n%g\n",  (1.0 - .2 - .2 - .2 - .2 - .2),
                       (1.0 - (.2 + .2 + .2 + .2 + .2))  );
   return 0;
}

5.55112e-017
0

Remarque hors-sujet (mais importante)

• Il y a pas mal de cas où ce qui change d'une expression à l'autre est le signe du zéro obtenu. Par exemple, x-y et -(y-x) ne sont pas équivalents car si les deux nombres sont égaux, alors la première expression donne +0.0 mais la deuxième -0.0.
• Il y a aussi des cas où les expressions ne sont pas équivalentes si on considère les valeurs "spéciales" des nombres, à savoir zéro, l'infini ou NaN :
  • x-x ne vaut pas 0.0 si x vaut l'infini ou NaN ;
  • 0*x ne vaut pas 0.0 si x vaut -0.0, l'infini ou NaN ;
  • x/x ne vaut pas 1.0 si x vaut zéro, l'infini ou NaN ;
  • x==x est faux si x vaut NaN ;
  • x!=x est vrai si x vaut NaN.

L'environnement des flottants : <fenv.h>

L'infâme traîtrise de l'opérateur ==

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12

int main(void) {
   float f=0;
   int i;
   
   for(i=0; i<100; ++i)
      f+= .01;
   
   if(f==1)   printf("f==1\n");
   else       printf("f!=1,  f==%f", f);
   
   return 0;
}

L'écart absolu & l'écart relatif

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12

#include <math.h>   // pour la fonction fabs, renvoyant la valeur absolue d'un flottant

#define  EPSILON  1e-8

/* ici pour des doubles, mais on fait la même chose pour des floats */
int doublesAreEqual(double a, double b) {
   return fabs(a-b) <= EPSILON;
/* La fonction fabs renvoie la valeur absolue (c'est-à-dire positive) du nombre
   de type double passé en argument ; ses équivalents (C99) sont fabsf pour les
   float, et fabsl pour les long double (pensez donc à adapter votre code en
   conséquence pour comparer les autres types flottants). */
}

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13

#include <math.h>

#define  EPSILON  1e-8

int doublesAreEqual(double a, double b) {
   if(a==b)   return 1;   // si a et b valent zéro  (explications plus bas)
   
   double absError= fabs(a-b);   // écart absolu
   a= fabs(a);   // on ne garde que les valeurs absolues pour la suite …
   b= fabs(b);   // … pour pouvoir calculer le nombre le plus grand en valeur absolue
   
   return  ( absError / (a>b? a:b) )  <=  EPSILON;
}

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13

#include <math.h>

#define  EPSILON  1e-8

int doublesAreEqual(double a, double b, double maxAbs, double maxRel) {
   if(a==b)   return 1;
   
   double absError= fabs(a-b);
   a= fabs(a);
   b= fabs(b);
   
   return   absError <= maxAbs   ||   ( absError / (a>b? a:b) )  <=  maxRel;
}

La représentation en mémoire convertie en entier

Valeur du float                       représentation en mémoire
                     binaire                                      hexadécimal    en base 10
+1.9999998           0   0111111 1   1111111 11111111 11111110    3F FF FF FE    1073741822
+1.9999999           0   0111111 1   1111111 11111111 11111111    3F FF FF FF    1073741823

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24

#include <stdlib.h>   // pour la fonction abs, renvoyant la valeur absolue d'un entier
#include <stdint.h>   /* header du C99, qui fournit des types entiers de taille fixe :
                         —  int32_t,  int64_t : entier signé de 32 bits ou de 64 bits ;
                         — uint32_t, uint64_t : entier non-signé de 32 bits ou de 64 bits. */

#define  INTREPOFFLOAT(f)   ( *(int32_t*)&(f) )   // représentation entière d'un float de 32 bits
#define  INTREPOFDOUBLE(d)  ( *(int64_t*)&(d) )   // représentation entière d'un double de 64 bits

#define  MAXULPS  5

/* ici pour des floats, mais on fait exactement pareil pour des doubles */
int floatsAreEqual(float a, float b) {
   if(a==b)   return 1;
   
   return abs( INTREPOFFLOAT(a) - INTREPOFFLOAT(b) )  <=  MAXULPS;
/* attention à la fonction abs, qui prend un int en argument ; un int fait 16
   ou 32 bits : il peut donc être trop petit pour contenir les 32 bits de la
   représentation entière d'un float, et sera de toutes façons insuffisant
   pour les 64 bits de celle d'un double. Voyez la fonction labs qui prend un
   long int, ou llabs (C99) qui prend un long long int, ou mieux, écrivez vos
   propres fonctions de valeur absolue, pour 32 et 64 bits (avec les types de
   <stdint.h>) ; ainsi, vous n'aurez plus de problèmes de taille des types
   pouvant varier. Ici, je garde les fonctions standards par souci de clarté. */
}

Valeur du float              représentation en mémoire
                         hexadécimal    en base 10 selon le complément à 2
+4.2038954 e-45          00 00 00 03     3
+2.8025969 e-45          00 00 00 02     2
+1.4012985 e-45          00 00 00 01     1
+0.0000000               00 00 00 00     0
-0.0000000               80 00 00 00    -2147483648
-1.4012985 e-45          80 00 00 01    -2147483647
-2.8025969 e-45          80 00 00 02    -2147483646
-4.2038954 e-45          80 00 00 03    -2147483645

Valeur du float              représentation "transformée"
                         hexadécimal    en base 10 selon le complément à 2
+4.2038954 e-45          00 00 00 03     3
+2.8025969 e-45          00 00 00 02     2
+1.4012985 e-45          00 00 00 01     1
+0.0000000               00 00 00 00     0
-0.0000000               00 00 00 00     0   *    (* = a été transformé)
-1.4012985 e-45          FF FF FF FF    -1   *
-2.8025969 e-45          FF FF FF FE    -2   *
-4.2038954 e-45          FF FF FF FD    -3   *

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23

#include <stdlib.h>
#include <stdint.h>

#define  INTREPOFFLOAT(f)   ( *(int32_t*)&(f) )
#define  INTREPOFDOUBLE(d)  ( *(int64_t*)&(d) )

#define  MAXULPS  5

int floatsAreEqual(float a, float b) {
   int32_t aInt= INTREPOFFLOAT(a);   // représentations entières
   int32_t bInt= INTREPOFFLOAT(b);
   
   if(aInt<0)   aInt= 0x80000000 - aInt;   // ou 0x8000000000000000 pour des doubles
   if(bInt<0)   bInt= 0x80000000 - bInt;
/* NOTE: on teste (aInt<0) et non (a<0). En effet, si a==-0.0 (zéro négatif),
   alors le test (a<0) renverrait faux, et on garderait la représentation de
   -0.0, à savoir 0x80000000. La règle du complément à 2 garde l'avantage du
   bit de signe : si le bit de poids fort est à 1, alors le nombre entier est
   négatif, et réciproquement ; on peut donc utiliser le test sur l'entier et
   non sur le flottant pour savoir s'il faut "transformer" la représentation. */
   
   return abs( aInt - bInt )  <=  MAXULPS;
}

• les infinis : comme vu précédemment, les infinis sont adjacents aux plus grands nombres en valeur absolue. Notre fonction pourrait par exemple nous dire qu'un nombre positif très très très grand et + sont égaux, alors que ce n'est pas vrai (ne discutez pas, d'un point de vue mathématique c'est faux ) ;
• les NaN : de même, certains NaN pourraient être comparés comme égaux à l'infini ou à un nombre extrêmement grand en valeur absolue, voire à un autre NaN. Or, normalement, n'importe quelle comparaison avec un NaN (hormis !=) devrait valoir faux.

Code complet

Header

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51

#ifndef INCLUDE_CMPFLOATS_H
#define INCLUDE_CMPFLOATS_H


#include <stdint.h>
#include <math.h>   /* pour les macros de test des nombres
                       flottants ( isnan() et isinf() ) */


/* représentation entière du nombre en virgule flottante */
#define  INTREPOFFLOAT(f)   ( *(int32_t*)&(f) )
#define  INTREPOFDOUBLE(d)  ( *(int64_t*)&(d) )

/* marge maximale séparant deux nombres flottants considérés comme égaux,
   en termes d'ULPs ("Unit of Least Precision") */
#define  MAXULPSFLOAT   5
#define  MAXULPSDOUBLE  5


/* renvoie vrai (1) si les 2 nombres sont égaux, faux (0) sinon */
int floatsAreEqual(float a, float b);
int doublesAreEqual(double a, double b);


/* renvoie  -1 si a>b,  0 si a==b,  1 si a<b,  ou -2 si a ou b est NaN */
int cmpFloats(float a, float b);
int cmpDoubles(double a, double b);

/* macros booléennes (à utiliser dans des tests simples) */
#define  CMPFLOATS_EQUAL(a,b)     (cmpFloats((a),(b))==0)      // =>  a==b
#define  CMPFLOATS_UNEQUAL(a,b)   (cmpFloats((a),(b))!=0)      // =>  a!=b
#define  CMPFLOATS_GT(a,b)        (cmpFloats((a),(b))==-1)     // =>  a>b
#define  CMPFLOATS_LT(a,b)        (cmpFloats((a),(b))==1)      // =>  a<b
//#define  CMPFLOATS_GTEQUAL(a,b)   (cmpFloats((a),(b))!=1)      // =>  a>=b
//#define  CMPFLOATS_LTEQUAL(a,b)   (cmpFloats((a),(b))!=-1)     // =>  a<=b
#define  CMPFLOATS_GTEQUAL(a,b)   ((cmpFloats((a),(b))-1)&2)   // =>  a>=b
#define  CMPFLOATS_LTEQUAL(a,b)   (cmpFloats((a),(b))>=0)      // =>  a<=b
#define  CMPFLOATS_NAN(a,b)       (cmpFloats((a),(b))==-2)     // =>  a==NaN || b==NaN

#define  CMPDOUBLES_EQUAL(a,b)     (cmpDoubles((a),(b))==0)    // =>  a==b
#define  CMPDOUBLES_UNEQUAL(a,b)   (cmpDoubles((a),(b))!=0)    // =>  a!=b
#define  CMPDOUBLES_GT(a,b)        (cmpDoubles((a),(b))==-1)   // =>  a>b
#define  CMPDOUBLES_LT(a,b)        (cmpDoubles((a),(b))==1)    // =>  a<b
//#define  CMPDOUBLES_GTEQUAL(a,b)   (cmpDoubles((a),(b))!=1)    // =>  a>=b
//#define  CMPDOUBLES_LTEQUAL(a,b)   (cmpDoubles((a),(b))!=-1)   // =>  a<=b
#define  CMPDOUBLES_GTEQUAL(a,b)   ((cmpDoubles((a),(b))-1)&2) // =>  a>=b
#define  CMPDOUBLES_LTEQUAL(a,b)   (cmpDoubles((a),(b))>=0)    // =>  a<=b
#define  CMPDOUBLES_NAN(a,b)       (cmpDoubles((a),(b))==-2)   // =>  a==NaN || b==NaN


#endif  //INCLUDE_CMPFLOATS_H

1
2
3
4
5

instructions1;
if(CMPFLOATS_GT(a,b)) { // a>b
   instructions2;
}
instructions3;

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12

instructions1;
if(CMPFLOATS_GT(a,b)) { // a>b
   instructions2;
}
else if(CMPFLOATS_LT(a,b)) { // a<b
   /*  /!\  on appelle la fonction cmpFloats 2 fois  */
   instructions2b;
}
else // a==b
   instructions2t;   /*  /!\  ce code est aussi exécuté dans le cas de NaN !  */
}
instructions3;

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14

instructions1;
switch(cmpFloats(a,b)) {
 case -1: // a>b
   instructions2;
   break;
 case  0: // a==b
   instructions2t;
   break;
 case  1: // a<b
   instructions2b;
   break;
 default: break; // NaN
}
instructions3;

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12

instructions1;
switch(cmpFloats(a,b)) {
 case -1: // a>b
 case  0: // a==b
   instructions2;   // ce code est donc exécuté si a>=b
   break;
 case  1: // a<b
   instructions2b;
   break;
 default: break; // NaN
}
instructions3;

Fichier source

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80

#include "cmpfloats.h"




/* valeur absolue d'un entier de 32 ou 64 bits */
uint32_t abs32(int32_t x) {   return x<0? -x : x;   }
uint64_t abs64(int64_t x) {   return x<0? -x : x;   }




int floatsAreEqual(float a, float b) {
   
   /* vérification pour NaN : si l'un des deux nombres est NaN, alors on
      retourne toujours faux */
   if(isnan(a) || isnan(b))
      return 0;
   
   /* vérification pour les infinis : si l'un des deux nombres est infini,
      alors on ne retourne vrai que si les deux nombres sont strictement
      égaux (tous les deux +inf ou -inf) */
   if(isinf(a) || isinf(b))
      return a==b;
   
   int32_t aInt= INTREPOFFLOAT(a);
   int32_t bInt= INTREPOFFLOAT(b);
   if(aInt<0)   aInt= 0x80000000 - aInt;
   if(bInt<0)   bInt= 0x80000000 - bInt;
   return abs32( aInt - bInt )  <=  MAXULPSFLOAT;
}

int doublesAreEqual(double a, double b) {
   if(isnan(a) || isnan(b))
      return 0;
   
   if(isinf(a) || isinf(b))
      return a==b;
   
   int64_t aInt= INTREPOFDOUBLE(a);
   int64_t bInt= INTREPOFDOUBLE(b);
   if(aInt<0)   aInt= 0x8000000000000000LL - aInt;
   if(bInt<0)   bInt= 0x8000000000000000LL - bInt;
   return abs64( aInt - bInt )  <=  MAXULPSDOUBLE;
}




int cmpFloats(float a, float b) {
   if(isnan(a) || isnan(b))
      return -2;   // -2 si on a au moins un NaN
   
   if(isinf(a) || isinf(b))
      return (a<b)? -1 : (a>b)? 1 : 0;   // gestion des infinis similaire à ci-dessous
   
   int32_t aInt= INTREPOFFLOAT(a);
   int32_t bInt= INTREPOFFLOAT(b);
   if(aInt<0)   aInt= 0x80000000 - aInt;
   if(bInt<0)   bInt= 0x80000000 - bInt;
   return (abs32(aInt-bInt) <= MAXULPSFLOAT)?  0   // 0 si les nombres sont égaux
                               : (aInt<bInt)?  1   // 1 si a<b
                               :              -1;  // -1 si a>b
}

int cmpDoubles(double a, double b) {
   if(isnan(a) || isnan(b))
      return -2;
   
   if(isinf(a) || isinf(b))
      return (a<b)? -1 : (a>b)? 1 : 0;
   
   int64_t aInt= INTREPOFDOUBLE(a);
   int64_t bInt= INTREPOFDOUBLE(b);
   if(aInt<0)   aInt= 0x8000000000000000LL - aInt;
   if(bInt<0)   bInt= 0x8000000000000000LL - bInt;
   return (abs64(aInt-bInt) <= MAXULPSDOUBLE)?  0
                                : (aInt<bInt)?  1
                                :              -1;
}

IEEE 754 et la norme C

En pratique : savoir si l'implémentation utilise IEEE 754

• dynamiquement (c'est-à-dire au moment de l'exécution de votre programme) : Vous pouvez par exemple déclarer un certain nombre de variables de type à virgule flottante, puis vérifier que leur représentation mémoire correspond à celle attendue en se basant sur IEEE 754. Il faudrait effectuer cette série de tests pour les différents types de nombres (zéros, dénormalisés, normalisés, infinis, NaNs).
  Un tel code ne serait pas infaillible (on peut très bien imaginer des formats ressemblant à IEEE 754, qui passeraient avec succès tous les tests) mais il permet d'éliminer certains formats.
  
  L'inconvénient est que du code inutile est intégré à l'exécutable, et que l'on perd du temps à chaque exécution du programme lorsqu'on effectue ces vérifications. À éviter en pratique, donc.
  D'un point de vue technique, cela reste cependant un bon exercice.
   
• statiquement (c'est-à-dire lors de la compilation) : Pour cela, il faut vous appuyer sur votre ami le préprocesseur. Sans détailler, vous pouvez tester la valeur des macros définies dans le header <float.h> du C99 ; celles-ci caractérisent l'implémentation des nombres flottants : précision, exposants minimum et maximum, valeurs minimales et maximales… Tout cela en 3 variantes pour chacun des 3 types à virgule flottante du C.
  
  Cette approche présente l'avantage de ne garder que le code nécessaire lors de la compilation et de ne pas faire cette vérification à l'exécution.
  Mais avouez que c'est fichtrement moins rigolo.