Strange results with floating-point comparison

Question

I have this simple test:

double h;
...
// code that assigns h its initial value, used below
...
if ((h>0) && (h<1)){
 //branch 1 -some computations
}
else{
 //branch 2- no computations
}

I listed my values as I got some really strange results and for example if: h=1 then the first branch is reached and I do not understand why since if h=1 I want branch2 to be computed.
Am I getting confused by something so obvious?

---

Edit:

This is how I compute and then use h:

double* QSweep::findIntersection(edge_t edge1,edge_t edge2) {  
point_t p1=myPoints_[edge1[0]];
point_t p2=myPoints_[edge1[1]];
point_t p3=myPoints_[edge2[0]];
point_t p4=myPoints_[edge2[1]];

double xD1,yD1,xD2,yD2,xD3,yD3,xP,yP,h,denom;
double* pt=new double[3];

// calculate differences  
xD1=p2[0]-p1[0];  
xD2=p4[0]-p3[0];  
yD1=p2[1]-p1[1];  
yD2=p4[1]-p3[1];  
xD3=p1[0]-p3[0];  
yD3=p1[1]-p3[1];    

xP=-yD1;
yP=xD1;
denom=xD2*(-yD1)+yD2*xD1;
if (denom==0) {
 return NULL;
}
else{
h=(xD3*(-yD1)+yD3*xD1)/denom;
}
std::cout<<"h is"<<h<<endl;
if (h < 1) std::cout<<"no"<<endl;
else std::cout<<"yes"<<endl;
if (h==1) {
 return NULL;
}
else{
if ((h>0)&&(h<1)){
 pt[0]=p3[0]+xD2*h;  
 pt[1]=p3[1]+yD2*h;
 pt[2]=0.00;
}
else{
 return NULL;
}
}


return pt;

}

---

Edit:

Okay, so it is clear how I should reformulate the condition.

From:

double h;
if (h==1){
   //computations here
}

To:

double h;
if (abs(h-1)<tolerance){
  //computations here
}

When I use double numbers.

But how do I reformulate this?

double h;
if (h<1){
   //computations here
}

Answer 1

Since h is a double, it may have been close enough to 1 to print as 1, but it is actually a bit less than 1 so the comparison succeeds. Floating-point numbers do that a lot.

Answer 2

It might be a precision issue. H might not be exactly 1, but something very near to it. Could you post some more information on what you're trying to do, for instance, where does the value of H come from, and where does it go?

Answer 3

if you have to use floats in checks, round them off and store it in for example a integer. 1f could be 1.0000000000000000000000000000000001 or 0.999999999999999999999999999999999999

Answer 4

Check the actual value of h by printing it out with maximum precision. You'll probably find that it is actually slightly less than 1.0.

I ran the following code as a test

#include <iostream>

int main()
{
    double h = 1.0;
    if((h>0) && (h<1))
    {
        std::cout << "first branch" << std::endl;
    }
    else
    {
        std::cout << "second branch" << std::endl;
    }
}

and the output was "first branch" (using g++ 4.3.2 on Ubuntu 8.10), but Indeera mentioned in a comment that the same code running on Windows XP compiled with VS2005 gives the output "second branch" (thanks, Indeera).

You might change your code to compare the differences between h and 0.0 and h and 1.0 to some small delta value.

double allowedDelta = 0.000001;

if( ((h - 0.0) > allowedDelta) && ((1.0 - h) > allowedDelta) )
... // h is between 0.000001 and 0.9999990

Note that "(h - 0.0)" can be replaced with "h" in this special case. I'm leaving it the way it is for illustrative value.

Also note that if you were only making one comparison you'd need to compare the delta to the absolute value of the difference between h and some constant. Since you're checking a range here, the two comparisons ANDed together make a special case where you can bypass the use of abs. If h is a negative value or some positive value greater than 1.0 it will be out of range and fail one of the two tests above.

Answer 5

It could have to do with the fact that doubles in C/C++ are 64bit, but computation could be done at higher precision (the floating point registers of your cpu are wider (96 bit), so that not even a statement like cos(x)==cos(x) could be true.

Reference: http://www.parashift.com/c++-faq-lite/newbie.html#faq-29.18

Answer 6

The reason is that floating point numbers are not a real representation of the number that you save in the variable. (Opposed to BCD [binary coded decimals])

You can see the definition here: Floating Point Wikipedia

So the problem is that certain numbers are not expressable with a given set of bits. (You could if you could add bits infinitely) The trick is that in most cases the difference between the saved number and the intended number are pretty small. In practice, you have some corner cases where this may lead to problems. This is for example the reason why you should not build a financial software and use floating point numbers for money calculations. You can easily have differences that are noticeable, which your tax office won't like.

So to compare floating point numbers you should always apply some kind of threshhold that is ok for your application. Something like:

if(a==b)

becomes

if(abs(a-b)<threshold)

Edit: As David Mentioned in his comment, you will still get into trouble with things like pi, 1./3., ... But you can at least store numbers without loss in precision that you put into the system. As computers have limited memory you can always construct corner cases where you can't rely on exact representations...

Just saw your edit of the text, so here is the next edit:

if(a<1)

Is somehow harder, because you don't know whether it's just the number representation that is wrong or if it's just a real value that is close to 1. It really depends on the requirements for your algorithm. If a small error is ok, then do:

if(a < 1-threshold)

If it's not ok, then you have to use another variable type that is not suffering from the problem.

Answer 7

[Removed, use comments or edit your question]

Answer 8

Always allow for rounding errors when comparing floating point values. Rather than testing for equality, something like this is usually what you want:

if (abs(x - y) < epsilon) ...

where epsilon is some suitably small value, like 0.00001.

There are several reasons why floating point math isn't accurate. First, not all values can be represented exactly (0.1 for example, can not be represented in binary floating point numbers, in the same way that 1/3 can't be represented in decimals)

Another is that floating point only uses a fixed number of significant digits (which "float" relative to the decimal point, hence the name). So on large numbers, small fractions effectively get truncated away. Never assume that floating point computations return an accurate result.

Answer 9

[Removed, use comments or edit your question]

Answer 10

Short story: Your tests are incorrect because floating point numbers do not behave like you probably expect them to. Particularly things like "denom == 0" are problematic.

Sun has been nice enough to provide this paper online:

What Every Computer Scientist Should Know About Floating Point Arithmetic

which is exactly as advertised. Depending on your background it will be an easy read or a big of work, but it really is worth the time for any programmer using floats.

ADDED COMMENT: I'm not suggesting that every programmer will easily understand everything in that paper. Reading over it though will at the very least give a better idea of what floats actually are, what the issues are, and some hints on how to handle things properly.

If you want to do a lot of numerical work properly, you're going to have to read up on the various techniques, but that would be a textbook (or several) worth of material. The comments here have already pointed out some of the basics, and linked to more

Answer 11

Okay, you've posted the code. You are calculating h by a series of arithmetic operations from what looks like fairly arbitrary numbers. This means you're going to get a very close approximation to the ideal value of h, but not quite the right one.

This means that you need to do approximate comparisons. Testing (h == 1.0) will succeed only by accident; try (fabs(h - 1.0) < 1e-10) or something like that (using a const double instead of a magic number for tolerance). Make the appropriate changes for other comparisons.

Answer 12

You may be very interested in the Numerical Robustness for Geometric Calculations (aka "EPSILON is NOT 0.00001!") talk given at GDC 2005.