Bugs in the FAQ....

[fungus]

I reported this bug years ago but it's still not fixed.

On this page: http://c-faq.com/fp/fpequal.html

it says to use the following code:

if (fabs(a - b) <= epsilon * fabs(a))

What happens when a is zero and b is (eg.) epsilon/2?

In this case a and b should compare as equal, but they don't...

 

[jacob navia]

I reported this bug years ago but it's still not fixed.

On this page: http://c-faq.com/fp/fpequal.html

it says to use the following code:

if (fabs(a - b) <= epsilon * fabs(a))

What happens when a is zero and b is (eg.) epsilon/2?

In this case a and b should compare as equal, but they don't...

If you read the next line it will say:

.... where epsilon is a value chosen to set the degree of ``closeness''
(and where you know that a will not be zero).

So a can't be zero

 

[fungus]

If you read the next line it will say:

... where epsilon is a value chosen to set the degree of ``closeness''
(and where you know that a will not be zero).

So a can't be zero

Ah, Ok. It didn't used to say that....

So is that the "fix" then?

 

[jacob navia]

Ah, Ok. It didn't used to say that....

So is that the "fix" then?

Well, it says that in the link you gave us.

Obviously it would be better if the code made no such an assumption like

if( a == 0 || fabs(a - b) <= epsilon * fabs(a))

but all the faqs have bugs and this one is not an exception. If you
reread the text it looks like the sentence within the parenthesis was
hastily added after someone discovered that...


well you know

 

[Willem]

fungus wrote:
) I reported this bug years ago but it's still not fixed.
)
) On this page: http://c-faq.com/fp/fpequal.html
)
) it says to use the following code:
)
) if (fabs(a - b) <= epsilon * fabs(a))
)
) What happens when a is zero and b is (eg.) epsilon/2?
)
) In this case a and b should compare as equal, but they don't...

Why should a and b compare as equal in that case ?
Consider the case where a = 0.001 and b = epsilon/2, for example.

However, a much more serious problem with the above code is that it
is not commutative. compare(a,b) is not always equal to compare(b,a).

To be commutative, shouldn't the code be something like:

if (fabs(a-b) <= epsilon * max(fabs(a), fabs(b)))

(Which, incidentally, solves the a=0 problem also)


SaSW, Willem
--
Disclaimer: I am in no way responsible for any of the statements
made in the above text. For all I know I might be
drugged or something..
No I'm not paranoid. You all think I'm paranoid, don't you !
#EOT

 

[fungus]

Obviously it would be better if the code made no such an assumption

I changed it to:

if (((a == 0) && (b<epsilon)) || (fabs(a - b) <= epsilon * fabs(a)))


There's still a _tiny_ bug because:

(epsilon * fabs(a)) != (epsilon * fabs(b))

But it's a lot better....

well you know

Yeah, I know...

 

[fungus]

Why should a and b compare as equal in that case ?

Because "epsilon" is pretty close to zero...?

 

[Keith Thompson]

Well, it says that in the link you gave us.

Obviously it would be better if the code made no such an assumption like

if( a == 0 || fabs(a - b) <= epsilon * fabs(a))

but all the faqs have bugs and this one is not an exception. If you
reread the text it looks like the sentence within the parenthesis was
hastily added after someone discovered that...

[...]

Yes, but that's not the right fix. It assumes that a and b are
"equal" whenever a==0.

 

[user923005]

I reported this bug years ago but it's still not fixed.

On this page:http://c-faq.com/fp/fpequal.html

it says to use the following code:

if (fabs(a - b) <= epsilon * fabs(a))

What happens when a is zero and b is (eg.) epsilon/2?

In this case a and b should compare as equal, but they don't...

It depends on whether you are trying to calculate relative or absolute
error.

This is one possible relative error comparison for floating point
types:

#include <float.h>
#include <math.h>

int double_compare (double d1, double d2)
{
if (d1 > d2)
if ((d1 - d2) < fabs (d1 * DBL_EPSILON))
return 0;
else
return 1;
if (d1 < d2)
if ((d2 - d1) < fabs (d2 * DBL_EPSILON))
return 0;
else
return -1;
return 0;
}

int float_compare (float d1, float d2)
{
if (d1 > d2)
if ((d1 - d2) < fabsf (d1 * FLT_EPSILON))
return 0;
else
return 1;
if (d1 < d2)
if ((d2 - d1) < fabsf (d2 * FLT_EPSILON))
return 0;
else
return -1;
return 0;
}

On the other hand, maybe you do not want relative error.
An example might be if you have a floating point function that is
given to a root search algorithm.
When the root has a y value of DBL_EPSILON you have found the zero and
you do not want to keep searching or refining. So the above method
would be incorrect for that application.
I don't think that there is a simple answer to the problem.