Implémenter assez_proches()

Quelques raccourcis :

Les équations sur ce site sont affichées avec MathJax.

Il arrive fréquemment que l'on souhaite implémenter une comparaison entre deux nombres à virgule flottante pour évaluer s'ils sont « égaux » ou non l'un à l'autre. Malheureusement, on confond souvent le domaine des réels , le domaine des rationnels et le domaine des nombres à virgule flottante (en C++, on parle de float, double et long double). En effet :

Exactitude et types

Concrètement, en programmation, c'est aux nombres à virgule flottante, aux entiers et aux rationnels que nous sommes typiquement confrontés. Les nombres à virgule flottante et les entiers ont droit à un support direct du substrat matériel, ce qui permet de les manipuler efficacement. Les entiers ont une représentation exacte, ce qui n'est en général pas le cas des nombres à virgule flottante comme nous venons de l'expliquer.

Empruntant un exemple à Bruce Dawson (article original), le programme suivant affichera probablement non.

#include <iostream>
int main() {
   using namespace std;
   float f = 1.1;
   if (f == 1.1)
      cout << "oui" << endl;
   else
      cout << "non" << endl;
}

La raison est simple : 1.1 est un littéral double, donc un nombre à virgule flottante à double précision, alors que f est un float, nombre à virgule flottante à simple précision. Puisque 1.1 n'a pas de représentation exacte sur un nombre à virgule flottante tel que représenté sur l'immense majorité (la totalité, probablement) des plateformes, alors le littéral 1.1 représente une approximation à environ une quinzaine de décimales près de , alors que f représente une approximation à environ sept ou huit décimales près de . Les deux ont une erreur, mais pas la même, donc l'évaluation de f == 1.1 nous donne faux.

En retour, les deux programmes ci-desous afficheraient oui : celui de droite parce que seuls des nombres à simple précision sont manipulés, et celui de gauche parce que seuls des nombres à double précision sont manipulés.

#include <iostream>
int main() {
   using namespace std;
   float f = 1.1f;
   if (f == 1.1f)
      cout << "oui" << endl;
   else
      cout << "non" << endl;
}
#include <iostream>
int main() {
   using namespace std;
   double f = 1.1;
   if (f == 1.1)
      cout << "oui" << endl;
   else
      cout << "non" << endl;
}

Enfin, examinons le programme suivant.

#include <iostream>
int main() {
   using namespace std;
   float f0, f1;
   cin >> f0 >> f1;
   float f = f0 - f1;
   cout << f0 << " - " << f1 << " == " << f << endl;
   if (f == 1.1f)
      cout << "oui" << endl;
   else
      cout << "non" << endl;
}

Ce programme lit deux float au clavier, calcule la différence f0-f1 et dépose le résultat dans f. Supposons qu'un individu entre au clavier 1.4 et 0.3 par exemple (la première ligne dans l'affichage ci-dessous), l'exécution donnera probablement ce qui suit :

1.4
0.3
1.4 - 0.3 == 1.1
non

Il peut apparaître surprenant que 1.1 sur des float soit différent du littéral 1.1f, aussi float.

Il faut cependant comprendre ici que f résulte d'un calcul sur des nombres approximatifs, et que ce calcul comprend fondamentalement une erreur, si petite soit-elle. Conséquemment, les deux 1.1 comparés avec == ne sont pas identiques si on les regarde dans le détail : ils diffèrent éventuellement d'une ou plusieurs décimales, et seuls les choix faits pour l'affichage des nombres à virgule flottante (par défaut, pour éviter la pollution à l'écran, les décimales les moins significatives ne sont pas affichées) nous cachent cette différence structurelle.

Cela dit, cette différence existe, et affecte nos calculs.

Premier constat : règle générale, on évitera de comparer des nombres à représentation inexacte avec == ou !=. Ceci est indépendant du langage de programmation, d'ailleurs. On sera aussi prudents avec les opérateurs relationnels que sont <, <=, > et >=.

Deuxième constat : en général, en comparant entre eux des nombres à virgule flottante, ce qui nous intéresse est surtout de savoir s'ils sont assez proches pour que la différence soit à l'intérieur de notre marge d'erreur, pas de savoir s'ils sont identiques (ceci n'étant en général pas pertinent, l'erreur étant intrinsèque à la représentation sous-jacente, quelle qu'elle soit).

Exemple simple : fonction assez_proches() naïve

Une implémentation simple d'une fonction assez_proches() sur deux float serait la suivante :

#include <cmath>
bool assez_proches(float f0, float f1, float epsilon = 0.00001f) {
   using std::abs;
   return abs(f0 - f1) <= epsilon;
}

Ce code n'est pas de qualité industrielle : il ne valide pas des cas où f0 serait le plus grand float possible (std::numeric_limits<float>::max()) et f1 serait non nul, par exemple, alors qu'un tel cas provoquerait un débordement de capacité.

La valeur est donnée par défaut à epsilon, pour éviter que le code client ne doive le suppléer en tout temps... Rien de nécessaire, mais c'est plus agréable à utiliser ainsi.

Il permet toutefois de comprendre le truc général : étant donné une certaine tolérance à l'erreur (donnée ici par le paramètre epsilon), on évalue si la valeur absolue de la différence entre les deux nombres est inférieure ou égale à cette tolérance. Dans ce cas, les deux nombres sont littéralement jugés assez proches l'un de l'autre.

Ce qui est agaçant avec cette implémentation est qu'il faut alors traiter les entiers et les nombres à virgule flottante différemment les uns des autres. Notez que j'ai utilisé les alternatives de C++ 17, qui permettent de localiser des variables à même un énoncé if, mais que ce n'est pas nécessaire à l'exemple (ce n'est qu'une question d'élégance) :

// ...
#include <iostream>
int main() {
   using namespace std;
   if (int a, b; cin >> a >> b)
      if (a == b)
         cout << a << " == " << b << endl;
      else
         cout << a << " != " << b << endl;
   if (float f0, f1; cin >> f0 >> f1)
      if (assez_proches(f0,f1))
         cout << a << " ~== " << b << endl;
      else
         cout << a << " ~!= " << b << endl;
}

Ceci peut poser problème dans du code générique. Par exemple :

//
// ceci est Ok pour un type T dont les valeurs sont
// représentées de manière exacte, mais ne l'est pas
// pour un type T comme float ou double
//
template <class T>
   void afficher_si_pas_trop_moche(const T &val, const T &moche, std::ostream &os) {
      if (val != moche) // <-- pas cool
         os << val;
   }

Une solution est de déterminer des variantes d'assez_proches() pour des types tels que int. Par exemple :

bool assez_proches(float f0, float f1, float epsilon = 0.00001f) {
   using std::abs;
   return abs(f0 - f1) <= epsilon;
}
bool assez_proches(int i0, int i1) {
   return i0 == i1;
}
//
// ceci est Ok pour tout type T pour lequel existe
// une définition de assez_proches(T,T)
//
template <class T>
   void afficher_si_pas_trop_moche(const T &val, const T &moche, std::ostream &os) {
      if (!assez_proches(val, moche))
         os << val;
   }

Dans la mesure où afficher_si_pas_trop_moche() ne tient pas compte d'un epsilon explicite, cette version fonctionnera à la fois pour float et pour int. Évidemment, couvrir explicitement tous les types entiers et tous les types à virgules flottante est fastidieux, alors nous souhaitons une solution plus générale.

Exemple plus général

Une implémentation un peu plus générale d'une fonction assez_proches(), tenant compte des types et se comportant différemment avec des nombres entiers et avec des nombres à virgule flottante, serait :

#include <limits>
#include <cmath>
class exacte {};
class approximative {};
template <bool>
   struct representation;
template <>
   struct representation<true> {
      using type = exacte;
   };
template <>
   struct representation<false> {
      using type = approximative;
   };
template <class T>
   bool assez_proches(const T &a, const T &b, approximative) {
      using std::abs;
      return abs(a-b) <= std::pow(T{10}, -5);
   }
template <class T>
   bool assez_proches(const T &a, const T &b, exacte) {
      return a == b;
   }
template <class T>
   bool assez_proches(const T &a, const T &b) {
      using std::numeric_limits;
      return assez_proches(
         a, b, typename representation<numeric_limits<T>::is_exact>::type{}
      );
   }

Ici, quelques techniques plus avancées ont été appliquées :

Un allègement syntaxique possible serait de créer un alias representation_t<bool> pour éviter le recours à typename ... ::type, comme suit. Cet allègement est conforme aux usages standardisés depuis C++ 14 :

#include <limits>
#include <cmath>
class exacte {};
class approximative {};
template <bool>
   struct representation;
template <>
   struct representation<true> {
      using type = exacte;
   };
template <>
   struct representation<false> {
      using type = approximative;
   };
template <bool V>
   using representation_t = typename representation<V>::type; // <-- ICI
template <class T>
   bool assez_proches(const T &a, const T &b, approximative) {
      using std::abs;
      return abs(a-b) <= std::pow(T{10}, -5);
   }
template <class T>
   bool assez_proches(const T &a, const T &b, exacte) {
      return a == b;
   }
template <class T>
   bool assez_proches(const T &a, const T &b) {
      using std::numeric_limits;
      return assez_proches(
         a, b, representation_t<numeric_limits<T>::is_exact>{} // <-- ICI
      );
   }

Ceci fonctionne bien mais a quelques défauts, incluant celui de ne pas permettre de spécifier la marge d'erreur pour les représentations inexactes – en particulier, dans ce cas, la même marge d'erreur s'appliquera à float, double et long double, et ce sera celle de la représentation la plus imprécise du lot, float, du fait que les marges d'erreurs appropriées pour les autres types ne pourront probablement pas être exprimées dans un float.

Diverses variantes

Diverses variantes sont possibles. En voici quelques-unes, en ordre grosso modo chronologique. En pratique, préférez la plus récente qui soit supportée par votre compilateur.

Variante avec traits (C++ 11)

Une variante, équivalente à bien des égards et utilisant des traits, serait (version C++ 11) :

#ifndef MATHS_EXT_H
#define MATHS_EXT_H
#include <type_traits>
#include <limits>
#include <cmath>
namespace maths_ext {
   //
   // ...
   // trucs divers dont on n'a pas besoin ici, omis
   // pour simplifier le portrait d'ensemble
   //
   template <class>
      struct numeric_limits_ext;
   template <>
      struct numeric_limits_ext<float> {
         static constexpr float seuil_assez_proches() noexcept {
            return 0.000001f;
         }
      };
   template <>
      struct numeric_limits_ext<double> {
         static constexpr double seuil_assez_proches() noexcept {
            return 0.0000000001;
         }
      };
   template <>
      struct numeric_limits_ext<long double> {
         static constexpr long double seuil_assez_proches() noexcept {
            return 0.000000000000001L;
         }
      };
   class exact_representation {};
   class float_representation {};
   template <class T>
      struct precision_traits {
         using representation_type = typename
            std::conditional<
               std::numeric_limits<T>::is_exact,
               exact_representation,
               float_representation
            >::type;
      };
   template <class T>
      bool assez_proches(const T &a, const T &b, exact_representation) {
         return a == b;
      }
   template <class T>
      bool assez_proches(const T &a, const T &b, float_representation) {
         return abs(a - b) <= numeric_limits_ext<T>::seuil_assez_proches();
      }
   template <class T>
      bool assez_proches(const T &a, const T &b) {
         return assez_proches(
            a, b, typename precision_traits<T>::representation_type{}
         );
      }
}
#endif

La technique est la même, mais j'ai ajouté un type numeric_limits_ext générique sur la base d'un type de nombre à virgule flottante, qui permet d'indique une marge d'erreur qui dépendra des types impliqués. J'ai aussi articulé la sélection de la catégorie de représentation sur la base d'un trait reposant sur une alternative statique, ce qui est plus général (ou, du moins, un peu moins artisanal).

Variante avec alias génériques (C++ 14)

On peut faire mieux encore, cela dit (code C++ 14), en utilisant des alias génériques et en maximisant le recours à constexpr :

#ifndef MATHS_EXT_H
#define MATHS_EXT_H
#include <type_traits>
#include <limits>
namespace maths_ext {
   //
   // ...
   // trucs divers dont on n'a pas besoin ici, omis
   // pour simplifier le portrait d'ensemble
   //
   template <class>
      struct numeric_limits_ext;
   template <>
      struct numeric_limits_ext<float> {
         static constexpr float seuil_assez_proches() noexcept {
            return 0.000001f;
         }
      };
   template <>
      struct numeric_limits_ext<double> {
         static constexpr double seuil_assez_proches() noexcept {
            return 0.0000000001;
         }
      };
   template <>
      struct numeric_limits_ext<long double> {
         static constexpr long double seuil_assez_proches() noexcept {
            return 0.000000000000001L;
         }
      };
   template <class T>
      constexpr T absolue(const T &val) {
         return val < 0? -val : val;
      }
   class exact_representation {};
   class float_representation {};
   template <class T>
      struct precision_traits {
         using type = std::conditional_t<std::numeric_limits<T>::is_exact, exact_representation, float_representation>;
      };
   template <class T>
      using precision_traits_t = typename precision_traits<T>::type;
   template <class T>
      constexpr bool assez_proches(const T &a, const T &b, exact_representation) {
         return a == b;
      }
   template <class T>
      constexpr bool assez_proches(const T &a, const T &b, float_representation) {
         return absolue(a - b) <= numeric_limits_ext<T>::seuil_assez_proches();
      }
   template <class T>
      constexpr bool assez_proches(const T &a, const T &b) {
         return assez_proches(a, b, precision_traits_t<T>{});
      }
}
#endif

Variante avec constante générique (C++ 17)

Depuis C++ 17, enfin, il est possible d'exprimer des variables génériques, ce qui permet de remplacer numeric_limits_ext<T>::seuil_assez_proches par une constante générique seuil_assez_proches<T>, dans le cas où le seuil serait le même pour toutes les implémentations mais où le type diffèrerait (en ce sens, nous aurions une implémentation légèrement différente de la précédente sur le plan sémantique) :

#ifndef MATHS_EXT_H
#define MATHS_EXT_H
#include <type_traits>
#include <limits>
namespace maths_ext {
   //
   // ...
   //
   template <class T>
      static constexpr T seuil_assez_proches = static_cast<T>(0.000001);
   template <class T>
      constexpr T absolue(const T &val) {
         return val < 0? -val : val;
      }
   class exact_representation {};
   class float_representation {};
   template <class T>
      struct precision_traits {
         using type = std::conditional_t<std::numeric_limits<T>::is_exact, exact_representation, float_representation>;
      };
   template <class T>
      using precision_traits_t = typename precision_traits<T>::type;
   template <class T>
      constexpr bool assez_proches(const T &a, const T &b, exact_representation) {
         return a == b;
      }
   template <class T>
      constexpr bool assez_proches(const T &a, const T &b, float_representation) {
         return absolue(a - b) <= seuil_assez_proches<T>;
      }
   template <class T>
      constexpr bool assez_proches(const T &a, const T &b) {
         return assez_proches(a, b, precision_traits_t<T>{});
      }
}
#endif

Variante avec enable_if

Une variante, un peu plus efficace en termes de temps de compilation, serait de remplacer la technique utilisée ici et reposant sur des classes vides à titre de types-étiquettes (pour distinguer les signatures de fonctions) par un recours à enable_if. Ce qui suit utilise C++ 17 par souci de concision et d'élégance, mais enable_if fonctionnait même avant C++ 11 :

#ifndef MATHS_EXT_H
#define MATHS_EXT_H
#include <type_traits>
#include <limits>
namespace maths_ext {
   //
   // ...
   //
   template <class T>
      static constexpr T seuil_assez_proches = static_cast<T>(0.000001);
   template <class T>
      constexpr T absolue(const T &val) {
         return val < 0? -val : val;
      }
   template <class T>
      constexpr std::enable_if_t<!std::is_floating_point_v<T>, bool> assez_proches(const T &a, const T &b) {
         return a == b;
      }
   template <class T>
      constexpr std::enable_if_t<std::is_floating_point_v<T>, bool> assez_proches(const T &a, const T &b) {
         return absolue(a - b) <= seuil_assez_proches<T>;
      }
}
#endif

Variante avec if constexpr

Avec l'avènement de C++ 17, nous avons désormais accès à if constexpr, ce qui permet d'alléger encore plus cette écriture :

#ifndef MATHS_EXT_H
#define MATHS_EXT_H
#include <type_traits>
#include <limits>
namespace maths_ext {
   //
   // ...
   //
   template <class T>
      static constexpr T seuil_assez_proches = static_cast<T>(0.000001);
   template <class T>
      constexpr T absolue(const T &val) {
         return val < 0? -val : val;
      }
   template <class T>
      constexpr bool assez_proches(const T &a, const T &b) {
         if constexpr(std::is_floating_point_v<T>)
            return absolue(a - b) <= seuil_assez_proches<T>;
         else
            return a == b;
      }
}
#endif

D'autres raffinements peuvent être envisagés encore.

Variante avec concepts

Depuis C++ 20, la programmation par concepts est une réalité, ce qui permet de réaliser des implémentations encore plus élégantes de fonctions telles que assez_proches(). Un concept permet d'exprimer des contraintes qui doivent être satisfaites par les types impliqués dans une fonction ou dans un type générique; programmer par concepts fait en quelque sorte de la programmation générique une programmation « normale », et offre les avantages de la programmation par interfaces sans en entraîner les inconvénients comme la perte de vitesse et le volet intrusif.

À cet effet, examinez ce qui suit :

#ifndef MATHS_EXT_H
#define MATHS_EXT_H
#include <type_traits>
#include <concepts>
namespace maths_ext {
   //
   // ...
   //
   template <class T>
      static constexpr T seuil_assez_proches = static_cast<T>(0.000001);
   template <class T> constexpr bool assez_proches(T a, T b) requires std::integral<T> {
      return a == b;
   }
   template <class T>
      constexpr T absolue(T val) { return val < 0? -val : val; }
   template <class T> constexpr bool assez_proches(T a, T b) requires std::floating_point<T> {
      return absolue(a - b) <= seuil_assez_proches<T>;
   }
#endif

L'expression requires indique la contrainte imposée au type: par exemple, requires std::integral<T> signifie qu'un fonction ne s'applique qu'aux types T qui satisfont le concept std::integral, ce dernier décrivant en termes formels ce qui est attendu d'un type entier. Ceci permet d'obtenir un effet semblable à ce qu'on pouvait réussir avec std::enable_if, mais de manière beaucoup plus élégante.

Notez que les expressions requires peuvent être complexes (p. ex. : requires C0<T> && C1<T> &&... pour les concepts C0, C1, ...), mais que nous n'avons pas besoin de cette complexité ici puisque dans les deux versions de notre fonction assez_proches(), nous n'imposons qu'une seule contrainte. À cet effet, examinez maintenant ce qui suit :

#ifndef MATHS_EXT_H
#define MATHS_EXT_H
#include <type_traits>
#include <concepts>
namespace maths_ext {
   //
   // ...
   //
   template <class T>
      static constexpr T seuil_assez_proches = static_cast<T>(0.000001);
   template <std::integral T> constexpr bool assez_proches(T a, T b) {
      return a == b;
   }
   template <class T>
      constexpr T absolue(T val) { return val < 0? -val : val; }
   template <std::floating_point T> constexpr bool assez_proches(T a, T b) {
      return absolue(a - b) <= seuil_assez_proches<T>;
   }
#endif

Ici, nous avons changé le terme général class (ou typename) des templates par un concept exprimant la contrainte imposée au type T. Ceci permet de la programmation générique par contrainte de manière simple et élégante.

Ici, pour les deux implémentations de la fonction assez_proches(), nous comparons respectivement deux nombres entiers ou deux nombres à virgule flottante de même type (donc deux int, deux short, deux double, etc. mais pas un float avec un double par exemple). Si nous souhaitions permettre la comparaison de deux variables de types distincts mais satisfaisant un même concept, nous aurions pu écrire :

#ifndef MATHS_EXT_H
#define MATHS_EXT_H
#include <type_traits>
#include <concepts>
namespace maths_ext {
   //
   // ...
   //
   template <class T>
      static constexpr T seuil_assez_proches = static_cast<T>(0.000001);
   template <std::integral T, std::integral U> constexpr bool assez_proches(T a, U b) {
      return a == b;
   }
   template <class T>
      constexpr T absolue(T val) { return val < 0? -val : val; }
   template <std::floating_point T, std::floating_point U> constexpr bool assez_proches(T a, U b) {
      return absolue(a - b) <= seuil_assez_proches<std::common_type_t<T,U>>;
   }
#endif

Notez le passage par std::common_type_t<T,U> pour choisir le type le plus adéquat pour le seuil de tolérance à l'erreur. Cette version permet donc de comparer un int avec un short, ou un float avec un double (dans quel cas, le seuil de tolérance à l'erreur sera de type double).

Si vous le souhaitez, il est possible de simplifier encore le tout. Examinez ce qui suit :

#ifndef MATHS_EXT_H
#define MATHS_EXT_H
#include <type_traits>
#include <concepts>
namespace maths_ext {
   //
   // ...
   //
   template <class T>
      static constexpr T seuil_assez_proches = static_cast<T>(0.000001);
   constexpr bool assez_proches(std::integral auto a, std::integral auto b) {
      return a == b;
   }
   template <class T>
      constexpr T absolue(T val) { return val < 0? -val : val; }
   constexpr bool assez_proches(floating_point auto a, floating_point auto b) {
      return absolue(a - b) <= seuil_assez_proches<std::common_type_t<decltype(a), decltype(b)>>;
   }
#endif

Ceci va au coeur de l'idée exprimée par notre algorithme.

Subtilités de l'arithmétique – il y a un bogue

Durant une présentation impromptue à CppNorth 2022 où j'ai utilisé cette technique pour résoudre une partie d'un probleme, on m'a signalé que cette solution peut mener à des bogues dans le cas d'un calcul impliquant un entier signé et un entier non-signé; particulièrement, assez_proches(-1,numeric_limits<unsigned int>::max()) retournera true dû à la promotion arithmétique qui convertira implicitement -1 en un immense unsigned int.

L'arithmétique, même des entiers, nous joue toutes et tous des tours un jour ou l'autre.

Il y a deux solutions à ce problème :

Une implémentation correcte serait donc :

#ifndef MATHS_EXT_H
#define MATHS_EXT_H
#include <type_traits>
#include <concepts>
#include <utility>
namespace maths_ext {
   //
   // ...
   //
   template <class T>
      static constexpr T seuil_assez_proches = static_cast<T>(0.000001);
   constexpr bool assez_proches(std::integral auto a, std::integral auto b) {
      return std::cmp_equal(a, b);
   }
   template <class T>
      constexpr T absolue(T val) { return val < 0? -val : val; }
   constexpr bool assez_proches(floating_point auto a, floating_point auto b) {
      return absolue(a - b) <= seuil_assez_proches<std::common_type_t<decltype(a), decltype(b)>>;
   }
#endif

Voilà!


Valid XHTML 1.0 Transitional

CSS Valide !