Alexis - 27 days ago 6
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
``````

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);
``````