28 Jun
2004
28 Jun
'04
8:06 a.m.
Hi all,
People who play around with floating point code (especially on x86)
quickly learn about the evils of comparing the equality of one floating
point value with another.
There are other related evils with floating point one of which I was
bitten by just recently and I thought I'd share it with y'all. If it
helps just one person from spending 20 odd hours chasing an elusive
bug like I did, I will have acheived something.
The evil I speak of is the difference between 32 and 64 bit floating
point representations (types float and double) and the x86 CPU's
internal 80 bit representation.
The most common trap is something like:
if (x == y)
do_something ();
where x and y are say double flotaing point numbers. Also lets just say
that the value of x is already in the CPU's FPU register as a result of
a previous operation and the other y, is not. What happens is that the
result of the previous operation can have a full 80 bits (part mantissa,
part exponent and a sign bit) of precision while y, loaded from memory
does not have this extra precision. The comparison therefore fails, even
though when printed out (or when compiler optimisation is switched off)
the two values are equal. This is the reason why the above if statement
is better written as:
if (fabs (x - y) < 1e-20)
do_something ();
The reason I am writing this email is that I was recently bitten by a
similar problem. I was keeping a running index into a table, and keeping
the integer part separate from the fractional part which was kept in a
double floating point:
double fractional = 0.0, increment = 0.1;
int integer = 0;
for (;;)
{
/* Bunch of other code. */
fractional += increment ;
integer += lrint (floor (fractional));
fractional -= floor (fractional);
}
The above code can produce very odd results for certain values of
increment. The problem in this case manifested itself in the
integer/fractional losing counts when compiled with gcc-3.4 while
the same code had worked perfectly with previous versions of the
compiler. The problem seems to be caused by the fact that the other
code in the loop was pushing at least some of the relevant values
out of the FPU stack into double floating point variables and that
when they were reloaded they had lost precision.
The fix in this case was this:
for (;;)
{
/* Bunch of other code. */
fractional += increment ;
rem = fmod (fractional, 1.0); /* floating point modulus */
integer += lrint (round (fractional - rem));
fractional = rem;
}
which is far more robust.
HTH,
Erik
--
+-----------------------------------------------------------+
Erik de Castro Lopo nospam(a)mega-nerd.com (Yes it's valid)
+-----------------------------------------------------------+
"The X-files is too optimistic. The truth is not out there."
-- Anthony Ord