Detect zero after std::remainder() call on doubles

Question

I have this code:

double a = remainder(92.66, 0.01);
double b = remainder(92.87, 0.01);

And it turns out that a = -5.33948e-15and b = -2.61423e-15

The answer here is clearly zero, and if I multiplied both numbers by 100 and did integer modulus it would be. But I'd like to be able to do this using doubles. The problem is that remainder is returning a number that is bigger than DBL_EPSILON - so what is the correct(best approximation) way to determine whether a or b is ~zero?

possible duplicate of [Is floating point math broken?](http://stackoverflow.com/questions/588004/is-floating-point-math-broken) — stefan, Jun 23 '14 at 12:46
Well maybe this "duplicate" isn't exactly the same, but come on: it's floating point maths all over again. There is no best way to compare floats. — stefan, Jun 23 '14 at 12:47
I understand how floating point arithmetic works, what I didn't understand is how to select a correct epsilon value in this case, which isn't in your duplicate. @Mark Ransom's answer was what I was looking for. Thanks. — Salgar, Jun 23 '14 at 13:20

score 2 · Answer 1 · answered Jun 23 '14 at 12:43

Comparing floating point numbers for equality is hard. Especially comparing them to zero. Here is an equality function copied directly from a useful blog post that you should read. Replace floats with doubles as you desire or create a template.

bool AlmostEqualRelativeAndAbs(float A, float B,
            float maxDiff, float maxRelDiff)
{
    // Check if the numbers are really close -- needed
    // when comparing numbers near zero.
    float diff = fabs(A - B);
    if (diff <= maxDiff)
        return true;

    A = fabs(A);
    B = fabs(B);
    float largest = (B > A) ? B : A;

    if (diff <= largest * maxRelDiff)
        return true;
    return false;
}

Now you need to figure out fitting values for the maxDiff and maxRelDiff parameters. They depend on how instable the calculation you perform is. Probably not much more than a few DBL_EPSILONs. Read more details in the linked article.

score 2 · Accepted Answer · answered Jun 23 '14 at 13:10

2

DBL_EPSILON is the smallest quantity that can be added to 1.0. If your number is larger than 1.0, the epsilon for that number will be proportionately larger as well. Simply multiply DBL_EPSILON by the number to get a threshold.

double a = remainder(92.66, 0.01);
if (a <= DBL_EPSILON * 92.66)
    a = 0.0;

answered Jun 23 '14 at 13:10

Mark Ransom

299,747
42
398
622

The first part is the part I didn't know. Thank you. – Salgar Jun 23 '14 at 13:19
The fundamental reason for the error is that none of your numbers can be represented by floats or doubles -- simple numbers like 0.1 are repeating decimals in binary. Since both 92.66 and 0.01 have some error the total error may be magnified by both and it takes careful analysis to figure out the bounds. If remainder is not correctly rounded (I don't know if correct rounding is guaranteed) then analysis gets even trickier. – Bruce Dawson Jun 23 '14 at 15:19

score 1 · Answer 3 · answered Jun 23 '14 at 13:01

Your problem is that your numbers are not exactly representable by the machine. You see the numbers 0.01, 92.66 and 92.87 in your code. The computer does not see them. It sees the closest numbers to 0.01, 92.66 and 92.87 that are represtable by the machine.

And when you ask it to compute the remainder, it computes the remainder of those representations and not of the original numbers you specified.

You must use one of many methods of comparing floating point to zero, depending on your application.

Detect zero after std::remainder() call on doubles

3 Answers3