#1 Wrong results for fcmp(1.0, (1.0+DBL_EPSILON)+DBL_EPSILON...

[Anonymous]

From k y n n a t p a n i x d o t c o m.

I'm not sure whether this bug is in the code or in the
documentation, so I will describe the issue in detail.

The README file gives this formula for choosing epsilon:

If the number of binary digits of error, e, is known,
then epsilon
can be calculated as follows:

epsilon = (pow(2, e) - 1) * FLT_EPSILON (for
floats)
epsilon = (pow(2, e) - 1) * DBL_EPSILON (for
doubles)

Suppose that I'm dealing with doubles and e == 1. I
interpret this to mean that I want to regard as "equal"
(i.e. "close enough") those doubles that differ only at
the last position; if they differ at a position other
than the last one they are to be considered different.

Since DBL_EPSILON is one ulp for doubles, 1.0 differs
from 1.0+DBL_EPSILON only at the last position (in
fact, this is how DBL_EPSILON is defined).
Furthermore, 1.0 differs from
(1.0+DBL_EPSILON)+DBL_EPSILON at the next-to-last
position. Moreover, according to the formula above,
for this case, epsilon should be DBL_EPSILON (since e
== 1).

Therefore, taking all this together, with this epislon

fcmp( 1., 1.+DBL_EPSILON, epsilon)

should be 0, and

fcmp( 1., (1.+DBL_EPSILON)+DBL_EPSILON, epsilon)

should be -1, but in fact, both expressions evaluate to
0. Unless I'm misunderstanding something, this looks
to me like a bug.

Here's some test code that illustrates what I just wrote:

#include <stdio.h>
#include <math.h>
#include <float.h>
#include "fcmp.h"

int main( void ) {
double x, y, z, w, epsilon;
x = 1.0;
y = x + DBL_EPSILON;
z = y + DBL_EPSILON;
w = z + DBL_EPSILON;
epsilon = DBL_EPSILON;

printf( "x %c= x\n", fcmp( x, x, epsilon ) ? '!' : '=' );
printf( "x %c= y\n", fcmp( x, y, epsilon ) ? '!' : '=' );
printf( "x %c= z\n", fcmp( x, z, epsilon ) ? '!' : '=' );
printf( "x %c= w\n", fcmp( x, w, epsilon ) ? '!' : '=' );

return 0;
}

If I compile and run this program, I get this output:

x == x
x == y
x == z
x != w

The puzzling result I mentioned earlier is the one on
the third line: "x == z".

Looking through the code for fcmp I see that the
results would agree with my expectations if I changed
the line in fcmp where delta is defined to

delta = ldexp(epsilon, exponent-1);

If I make this change, then the output of the test
program becomes

x == x
x == y
x != z
x != w

...which makes sense to me.

Alternatively, one could change the formula for
deriving epsilon to

epsilon = ((pow(2, e) - 1) * FLT_EPSILON)/2 (for
floats)
epsilon = ((pow(2, e) - 1) * DBL_EPSILON)/2 (for
doubles)