C++ Comparison of doubles using an epsilon

Question

I defined a comparison function for doubles that uses an epsilon, as suggested in: What is the most effective way for float and double comparison?. Two doubles differ if their absolute difference is less than epsilon.

If I have doubles d1 and d2, with d2 = d1+epsilon, d1 should not be equal to d2 since they differ by epsilon. Right?

It works with some values, but not all the time. It might be related to the architecture or to the compiler? What should I do to improve my comparison function?

#include <iomanip>  // For std::setprecision
#include <iostream>
#include "math.h"   // For fabs

const double epsilon = 0.001;

bool equal(const double &d1, const double &d2)
{
    return fabs(d1-d2) < epsilon;
}

int main()
{
    double d1 = 3.5;
    double d2 = d1+epsilon;

    // d2 = d1 + epsilon should not be equal to d1.
    std::cout << std::setprecision(16) << "d1 = " << d1 << " d2 = " << d2 << std::endl;
    bool eq = equal(d1,d2);
    if (eq)
    {
        std::cout << "They should not be equal, but they are!" << std::endl;
    }
    else
    {
        std::cout << "They are not equal. OK!" << std::endl;
    }
}

Output on my machine:

d1 = 3.5 d2 = 3.501

They should not be equal, but they are!

"If I have doubles d1 and d2, with d2 = d1+epsilon, d1 should not be equal to d2 since they differ by epsilon. Right?" No, exactly by the same reason. — Slava, Dec 22 '17 at 15:33
As a side-note the `std::setprecision` stream modifier doesn't modify the built-in type precision. You can't modify the built-in type precision. — Ron, Dec 22 '17 at 15:34
Don't you see the irony? You use epsilon comparison in order to combat the fact that floating point arithmetic is imprecise - then expect precise results for expressions involving `epsilon`. The exact mathematical value of `d1+epsilon` is not representable in `double`, so `d2` ends up being slightly smaller. — Igor Tandetnik, Dec 22 '17 at 15:35

score 2 · Answer 1 · answered Dec 22 '17 at 15:42

2

One of the most frequent floating-point caveats: if you put a number to the floating point, it doesn't mean it will be stored as exactly the same value. Instead, they would be slightly amended to create a binary structure that is most close to your intended value.

Ok, back to business! Try to add this to your source code:

std::cout << fabs(d1-d2) << std::endl;

And find the output as the following :)

0.0009999999999998899

You see? It's clearly less than your epsilon.

answered Dec 22 '17 at 15:42

Yury Schkatula

5,291
2
18
42

I was expecting no precision issue because of the "large" value of epsilon (0.001) compared to the system precision for doubles (~1e-16), clearly that was wishful thinking. – Montplaisir Dec 22 '17 at 16:16
This is regardless the system precision. Floating-point numbers are just limited amount of bits that encode a value. Obviously, the amount of floating-point values is non-countable and amount of bit combinations is countable. So, any regular floating point format (float, double and so on) is loss-driven. So, try to think like that: floating-points are kind of JPEGs (you get storage-based artifacts anyway) and ints are PNGs :) – Yury Schkatula Dec 22 '17 at 16:28

Slava · Answer 2 · 2017-12-22T15:58:57.593

If I have doubles d1 and d2, with d2 = d1+epsilon, d1 should not be equal to d2 since they differ by epsilon. Right?

Wrong. Your logic is broken. You've told that you should not compare doubles without epsilon but next step you are surprised that fabs(d1-d2) == epsilon is not true.

What should I do to improve my comparison function?

Nothing, d1 + epsilon is a border case and you cannot be sure if that number is equal to d1 or not. That is downside of floating point arithmetic.

score 1 · Answer 3 · answered Dec 22 '17 at 16:08

There are actually two different problems here. First, your comparison should be >=, rather than >, so it “works” for your test case where the difference is exactly epsilon.

The second (potential) problem is that you can’t really use a constant epsilon for these sorts of comparisons, in the general case. You need to select the epsilon based on the expected magnitude of the numbers involved, and the expected accumulated numerical error in the calculations.

Floating-point numbers are less precise the larger the magnitude of the value they’re storing is. If you repeat your test with a value that’s “large enough”, it’ll fail every time.

Rama · Answer 4 · 2017-12-22T18:21:26.773

Let's call pe precision error. It is infinitesimal(~2^−53) for IEEE 754 standard binary64

The maths are: [Pseudo code]

d1 = 3.5 + pe1
d2 = 3.501 + pe2

d2 - d1 = 0.001 + pe2 - pe1

epsilon = 0.001 + pe3

The result of the expression: [Pseudo code]

|d1-d2| < epsilon  //Since you know that d2>d1 you can replace |d1-d2| by d2-d1
0.001 + pe2 - pe1 < 0.001 + pe3
pe2 - pe1 < pe3

This border case can be true or false dependig on magnitude and sign of pe3, pe2 and pe1. But it is undetermined.

C++ Comparison of doubles using an epsilon

4 Answers4