Alexis Alexis - 2 months ago 10
C Question

Comparing doubles yields one less AND equal to the other

So I have the following code:

double which_min(const int i, const int j, const int nx,
const double step, const double local_cost,
double *D)
{
double tuple[3];

// DIAG, LEFT, UP
tuple[0] = (D[d2s(i-1, j-1, nx)] == NOT_VISITED) ? DBL_MAX : D[d2s(i-1, j-1, nx)] + step * local_cost;
tuple[1] = (D[d2s(i, j-1, nx)] == NOT_VISITED) ? DBL_MAX : D[d2s(i, j-1, nx)] + local_cost;
tuple[2] = (D[d2s(i-1, j, nx)] == NOT_VISITED) ? DBL_MAX : D[d2s(i-1, j, nx)] + local_cost;
/*
if (i == 83 && j == 124) printf("less? %d\n", tuple[1] < tuple[0]);
if (i == 83 && j == 124) printf("equal? %d\n", tuple[1] == tuple[0]);
if (i == 83 && j == 124) printf("greater? %d\n", tuple[1] > tuple[0]);
*/
int min = (tuple[0] <= tuple[1]) ? 0 : 1;
min = (tuple[min] <= tuple[2]) ? min : 2;

if (i == 83 && j == 124) printf("min = %d\n", min + 1);

return ((double) min + 1.0);
}


For the problematic case I'm working with, when
step = 1
,
i = 83
and
j = 124
, I'm getting

min = 2.


However, if I uncomment the other
if
statements, I get the correct value of
min
, but...

less? 1
equal? 1
greater? 0
min = 1


I know comparing doubles can be finicky, and I read this answer that mentions moving from 80-bit register to 64-bit memory, so I'm guessing the
printf
statements could be causing something like that. However,


  1. It still feels counter intuitive that
    tuple[1]
    is both less AND equal to
    tuple[0]
    , how is that possible?*

  2. Can I change this so that it gives the same result consistently?



And by the way, this only happens in a 32-bit system, the 64-bit version doesn't give me any problems.

EDIT: I guess that the first
printf
results in a loss of precision so that the
==
comparison evaluates to
1
afterwards, but still, the main question is whether I can make the result consistent.

EDIT2: Indeed, changing the order of
if
statements to

if (i == 83 && j == 124)
{
printf("equal? %d\n", tuple[1] == tuple[0]);
printf("less? %d\n", tuple[1] < tuple[0]);
printf("greater? %d\n", tuple[1] > tuple[0]);
}


Results in

equal? 0
less? 0
greater? 0
min = 1


And the values I get with
printf(%a)
are

tuple[0] = 0x1.2594a8056d275p+5
tuple[1] = 0x1.2594a8056d275p+5
tuple[2] = 0x1.fffffffffffffp+1023

Answer

Indeed the results stored in tuple were being read with more precision than expected, causing some of the comparisons to yield contradicting results.

Adding the -ffloat-store flag to the compiler remedies the situation, but as mentioned in the comments, it might not be ideal, and apparently it isn't portable.

Declaring tuple as volatile seems to do the trick, but it had to be done as

volatile double *tuple = (double *)malloc(3 * sizeof(double));

And then

free((double *)tuple);