CODEWITHSUNDEEP
c++floating-point
Guillaume Paris
Guillaume Paris

Reputation: 10539

c++ comparing two floating point values

I'm wondering what is the difference for comparing two double between this two manner :

double a1 = ...;
double a2 = ....;
  1. fabs(a1-a2) < epsilon
  2. (fabs(a1-a2)/a2) < epsilon

Is there a prefering way to do that ?

thanks

Upvotes: 5

Views: 16206

Answers (5)

Daniel La&#252;gt
Daniel La&#252;gt

Reputation: 1147

Epsilon varies according to the value inside the double range. If you want to use your own and fixed epsilon, I recommand to choose an epsilon for 1.

fabs(a1-a2) < epsilon

This is not perfect. If a1 and a2 are results from operations with small numbers, the epsilon would be small. If a1 and a2 are results from operations with big numbers, the epsilon would be big.

fabs((a1-a2)/a2) < epsilon

This is a bit better as you would like to scale epsilon by a2. However, there is a division by 0 if a2 equals 0.

fabs(a1-a2) < fabs(a2)*epsilon

This is a bit better. However, this is incorrect if a2 equals 0. fabs(a2)*epsilon would be equal to 0 and the comparison of two values a1=0.000000001 and a2=0 will always fail.

fabs(a1-a2) < max(fabs(a1), fabs(a2))*epsilon

This is a bit better. However, this is incorrect because epsilon is not proportional in continuous manner. With IEEE encoding, epsilon is proportional to a value in discrete manner on base 2.

fabs(a1-a2) < 2^((int)log2(max(fabs(a1), fabs(a2)))*epsilon

This looks correct for me when both a1 and a2 are not equal to 0.

The implementation of a generic method could be:

bool nearly_equal(double a1, double a2, double epsilon)
{
  if (a1 == 0 && a2 == 0)
    return true;

  return std::abs(a1 - a2) < epsilon * pow (2.0, static_cast<int> (std::log2(std::max(std::abs(a1), std::abs(a2)))));
}

Upvotes: 3

Daniel La&#252;gt
Daniel La&#252;gt

Reputation: 1147

Personally, I'm using std::nextafter for comparing two double. This use the smallest epsilon on a value (or a factor of the smallest epsilon). 

bool nearly_equal(double a, double b)
{
  return std::nextafter(a, std::numeric_limits<double>::lowest()) <= b
    && std::nextafter(a, std::numeric_limits<double>::max()) >= b;
}

bool nearly_equal(double a, double b, int factor /* a factor of epsilon */)
{
  double min_a = a - (a - std::nextafter(a, std::numeric_limits<double>::lowest())) * factor;
  double max_a = a + (std::nextafter(a, std::numeric_limits<double>::max()) - a) * factor;

  return min_a <= b && max_a >= b;
}

Upvotes: 4

Martin Stone
Martin Stone

Reputation: 13007

This article answers your question quite thoroughly, I think. You might want to skip ahead to the section "Epsilon comparisons".

Upvotes: 13

Stephen Canon
Stephen Canon

Reputation: 106167

This has already been pretty well addressed, but a few points are in order:

fabs(a1-a2) < epsilon

is comparing the absolute difference between a1 and a2 to the tolerance epsilon. This may be appropriate if you know the scaling a priori (for example, if a2 is actually a constant), but should generally be avoided if you don't know how big a1 and a2 are.

Your second option almost computes the relative difference, but has a bug; it should actually read:

fabs((a1-a2)/a2) < epsilon

(note that the division is inside the absolute value; otherwise, this condition is useless for negative a2). Relative error is more correct for most uses, because it more closely mirrors the way floating-point rounding actually happens, but there are situations in which it does not work and you need to use an absolute tolerance (usually this is because of catastrophic cancellation). You also will sometimes see relative error bounds written in this form:

fabs(a1-a2) < fabs(a2)*epsilon

which is often somewhat more efficient because it avoids a division.

Upvotes: 2

Pontus Gagge
Pontus Gagge

Reputation: 17258

The former only compares the absolute values, whereas the second compares the relative values. Say that epsilon is set to 0.1: this may be sufficient if a1 and a2 are largish. However, if both values are close to zero, the first way will consider most values equal.

It really depends on what kinds of values you are dealing with. Just be sure to consider the case a2==0 if you use the somewhat more mathematically reasonable second case.

Upvotes: 2

Related Questions