Comparing FloatsHow To Determine if Floating Quantities Are Close Enough Once a Tolerance Has Been Reached

• By Alberto Squassabia
• March 16, 2000

#ifndef _CMPFLT_HXX
#define _CMPFLT_HXX

#include <limits>
#include <float.h>
#include <stdlib.h>

//
// your favorite tolerances here
//

#define _FLT_TOL_MULT     10.0F
#define _DBL_TOL_MULT    100.0
#define _LDBL_TOL_MULT   1000.0L

// -------------------------------------
// inexact types tolerances start
// -------------------------------------

template <class nType> struct nTol {
   
   //
   // used only for exact types at compilation time; never run
   // except if floating point specialization is missing.  If so,
   // abort() is invoked.
   //
   typedef nType typedQty;

   static typedQty TolAtEPS() {
      if ( ! std::numeric_limits<nType>::is_exact ) abort();
         return 1;
   };
   static typedQty TolAtMIN() { 
      if ( ! std::numeric_limits<nType>::is_exact ) abort();
         return 1; 
   };
   static typedQty onePlusTol()  {
      if ( ! std::numeric_limits<nType>::is_exact ) abort();
         return 1; 
   };
   static typedQty oneMinusTol() { 
      if ( ! std::numeric_limits<nType>::is_exact ) abort();
         return 1; 
   };
   static typedQty zero()        { 
      if ( ! std::numeric_limits<nType>::is_exact ) abort();
         return 0; 
   };
   static typedQty one()         { 
      if ( ! std::numeric_limits<nType>::is_exact ) abort();
         return 1; 
   };
};

struct nTol<float> {
// -------------------------------------


   typedef float typedQty;

   static inline typedQty TolAtEPS()	{ return FLT_EPSILON *                   
	 _FLT_TOL_MULT; };
   static inline typedQty TolAtMIN()	{ return FLT_MIN *
	 _FLT_TOL_MULT; };
   static inline typedQty onePlusTol()	{ return 1.0F + FLT_EPSILON * 
	 _FLT_TOL_MULT; };
   static inline typedQty oneMinusTol()	{ return 1.0F - FLT_EPSILON * 
	 _FLT_TOL_MULT; };
   static inline typedQty zero()	{ return 0.0F; };
   static inline typedQty one()	{ return 1.0F; };
};

struct nTol<double> {
// -------------------------------------

   typedef double typedQty;

   static inline typedQty TolAtEPS()	{ return DBL_EPSILON * 
	 _DBL_TOL_MULT; };
   static inline typedQty TolAtMIN()	{ return DBL_MIN *
	 _DBL_TOL_MULT; };
   static inline typedQty onePlusTol()	{ return 1.0 + DBL_EPSILON * 
	 _DBL_TOL_MULT; };
   static inline typedQty oneMinusTol()	{ return 1.0 - DBL_EPSILON * 
	 _DBL_TOL_MULT; };
   static inline typedQty zero()	{ return 0.0; };
   static inline typedQty one()	{ return 1.0; };
};

struct nTol<long double> {
// -------------------------------------

   typedef long double typedQty;

   static inline typedQty TolAtEPS()	{ return LDBL_EPSILON * 
	 _LDBL_TOL_MULT; };
   static inline typedQty TolAtMIN()	{ return LDBL_MIN * 
	 _LDBL_TOL_MULT; };
   static inline typedQty onePlusTol()	{ return 1.0L + LDBL_EPSILON * 
	 _LDBL_TOL_MULT; };
   static inline typedQty oneMinusTol()	{ return 1.0L - LDBL_EPSILON * 
	 _LDBL_TOL_MULT; };
   static inline typedQty zero()       	{ return 0.0L; };
   static inline typedQty one()         	{ return 1.0L; };
};


// -------------------------------------
// inexact types tolerances end
// -------------------------------------


template <class nType> bool testNumEq (const nType& a, const nType& b) {
// -------------------------------------
// type-safe numerical equality
// -------------------------------------

   //
   // easy way out
   //

   if (a == b) return true;
   // check if exact (i.e. integer type); if so, equality failed
   // (just in case called with exact type, not a likely event)
   if ( std::numeric_limits<nType>::is_exact ) return false;
//
   // machinery for "close enough"
   //

   // type identification, storage
   typedef nTol<nType> typedQty;
   nType ratio, fabsa, fabsb;

   // test for closeness to 0
   if ( typedQty::zero() == a ) // check if |b| is "close enough" to 0 
      return (typedQty::TolAtEPS() > (typedQty::zero() > b ? -b : b));  
   if ( typedQty::zero() == b ) // check if |a| is "close enough" to 0 
      return (typedQty::TolAtEPS() > (typedQty::zero() > a ? -a : a));  

   // test for closeness to each other
   if (typedQty::TolAtMIN() > (fabsb=(typedQty::zero() > b ? -b : b))) { 
      // |b| is undistinguishable from 0
      if (typedQty::TolAtMIN() > (typedQty::zero() > a ? -a : a)) return true; 
      // . . . but |a| is distinguishable from 0
      return false;
   } else if (typedQty::TolAtMIN() < (fabsa=(typedQty::zero() > a ? -a : a))) {
      // |a| and |b| are distinguishable from 0, take ratio 
      // avoid overflow if |a| is very large and |b| is very small
      if ( (fabsb < typedQty::one()) 
        && (fabsa > fabsb * std::numeric_limits<nType>::max()) ) return false;
      // avoid underflow if |a| is very small and |b| is very large 
      if ( (fabsb > typedQty::one()) 
        && (fabsa < fabsb * std::numeric_limits<nType>::min()) ) return false;
      // this must be signed
      ratio = a/b;
   } else {
      // |b| is distinguishable from 0
      // |a| is undistinguishable from 0 
      return false;
   };

   // only if signed ratio is indistinguishable from +1 then equality true
   if ( (typedQty::onePlusTol() > ratio) 
         && (typedQty::oneMinusTol() < ratio)) return true;

   // they're different 
   return false;
};

#endif