When is a>a true ?

Question

Right, I think I really am living a dream. I have the following piece of code which I compile and run on an AIX machine:

AIX 3 5
PowerPC_POWER5 processor type
IBM XL C/C++ for AIX, V10.1
Version: 10.01.0000.0003


#include <stdio.h>
#include <math.h>

#define RADIAN(x) ((x) * acos(0.0) / 90.0)

double nearest_distance(double radius,double lon1, double lat1, double lon2, double lat2){
    double rlat1=RADIAN(lat1);
    double rlat2=RADIAN(lat2);
    double rlon1=lon1;
    double rlon2=lon2;
    double a=0,b=0,c=0;

    a = sin(rlat1)*sin(rlat2)+ cos(rlat1)*cos(rlat2)*cos(rlon2-rlon1);
    printf("%lf\n",a);
    if (a > 1) {
      printf("aaaaaaaaaaaaaaaa\n");
    }
    b = acos(a);
    c = radius * b;

    return radius*(acos(sin(rlat1)*sin(rlat2)+
        cos(rlat1)*cos(rlat2)*cos(rlon2-rlon1)));

}

int main(int argc, char** argv) {
  nearest_distance(6367.47,10,64,10,64);
  return 0;
}

Now, the value of 'a' after the calculation is reported as being '1'. And, on this AIX machine, it looks like 1 > 1 is true as my 'if' is entered !!! And my acos of what I think is '1' returns NanQ since 1 is bigger than 1. May I ask how that is even possible ? I do not know what to think anymore !

The code works just fine on other architectures where 'a' really takes the value of what I think is 1 and acos(a) is 0.

Answer 1

Print out the bits. You might just be getting fooled by some rounding error in the display of the floats as decimal real numbers.

Answer 2

1.000000000000000000001 is greater than 1. Are you sure you just aren't see enough decimal places? If that check is passing I'd wager thats your issue.

The usual solution is to use some form of epsilon to stop you worrying about rounding errors. ie if the double you have ought to be then try doing

if ( a > 1.00001f )

Its probably close enough to one so as not to cause you problems :)

Answer 3

The printf function, without a specified precision, will only show you the first 6 digits. So, try printing with a higher degree of precision... it is possible that a is slightly larger than 1, but only by a little. If you want to make things more robust, instead of (a>1), you can use (a-1)>epsilon for some value of epsilon.

Answer 4

If you do a comparison where result and expctedResult are float types:

if (result == expectedResult)

Then it is unlikely that the comparison will be true. If the comparison is true then it is probably unstable – tiny changes in the input values, compiler, or CPU may change the result and make the comparison be false.

Comparing with epsilon – absolute error

if (fabs(result - expectedResult) < 0.00001)

From Comparing floating point numbers

---

What Every Computer Scientist Should Know About Floating-Point Arithmetic