Floating point behavior

Floating point behavior

[eool]

Hi all,

I was reading the C++ FAQ at http://www.parashift.com/c++-faq-lite/ when I
encountered an interesting curiosity in the "newbie questions" section. The
FAQ warns against directly comparing floating point numbers and, in section
29.18, gives an example of code (see below) that produces behavior that
seems unnatural to many (most?) programmers. (Below code copied from
http://www.parashift.com/c++-faq-lite/newbie.html#faq-29.18)

 #include <cmath>
 
 void foo(double x, double y)
 {
   if (cos(x) != cos(y)) {
     std::cout << "Huh?!?\n";  ← you might end up here when x == y!!
   }
 }
 
 int main()
 {
   foo(1.0, 1.0);
   return 0;
 } 

Basically, the reason given in the FAQ is that the result of cos(x)
(intermediate, higher precision value) will be stored in a register and thus
truncated. Then the result of cos(y) is directly compared with the truncated
value located in a register, resulting in this strange (?) behavior. Now, my
question is, why isn't gcc able to tell that it is about to compare two
floating point numbers, one of which has been truncated, and then similarly
truncate the second value before the comparison? (In fact, since the return
type of cos() is double, shouldn't the intermediate value be truncated to
fit a double in any case?)

Second (this question might be even more naive than the previous one), if
you are never, ever, supposed to compare floating point numbers using ==; if
you are always supposed to check the difference of the numbers against an
epsilon value; then how come the compiler can't automatically do this? Is it
because analyzing the source code to track the accumulated floating point
errors (to appropriately adjust epsilon) might be prohibitively
difficult/slow? Or are there actually any situations where you'd want to
compare two floating point numbers directly with each other?

Thanks for taking the time to answer my "newbie questions" :)

/Erik
-- 
View this message in context: http://www.nabble.com/Floating-point-behavior-tp19304865p19304865.html
Sent from the gcc - Help mailing list archive at Nabble.com.

Re: Floating point behavior

[Andrew Haley]

eool wrote:

> I was reading the C++ FAQ at http://www.parashift.com/c++-faq-lite/ when I
> encountered an interesting curiosity in the "newbie questions" section. The
> FAQ warns against directly comparing floating point numbers and, in section
> 29.18, gives an example of code (see below) that produces behavior that
> seems unnatural to many (most?) programmers. (Below code copied from
> http://www.parashift.com/c++-faq-lite/newbie.html#faq-29.18)
> 
>  #include <cmath>
>  
>  void foo(double x, double y)
>  {
>    if (cos(x) != cos(y)) {
>      std::cout << "Huh?!?\n";  ← you might end up here when x == y!!
>    }
>  }
>  
>  int main()
>  {
>    foo(1.0, 1.0);
>    return 0;
>  } 
> 
> Basically, the reason given in the FAQ is that the result of cos(x)
> (intermediate, higher precision value) will be stored in a register and thus
> truncated. Then the result of cos(y) is directly compared with the truncated
> value located in a register, resulting in this strange (?) behavior. 

The results of Intel's builtin trig functions are long double, not double.
In the example above, the first result may be stored in memory as a double --
therefore truncated -- and compared with a direct result that is long
double.

> Now, my
> question is, why isn't gcc able to tell that it is about to compare two
> floating point numbers, one of which has been truncated, and then similarly
> truncate the second value before the comparison? (In fact, since the return
> type of cos() is double, shouldn't the intermediate value be truncated to
> fit a double in any case?)

I suppose it might be possible, but it would require considerable work.

> Second (this question might be even more naive than the previous one), if
> you are never, ever, supposed to compare floating point numbers using ==; if
> you are always supposed to check the difference of the numbers against an
> epsilon value; then how come the compiler can't automatically do this?

Three reasons:

1.  Anyone with a clue does this anyway.
2.  We compiler writers have no idea what an appropriate value of Epsilon
might be for your application.
3.  The standard doesn't allow us to change code in this way.

> Is it
> because analyzing the source code to track the accumulated floating point
> errors (to appropriately adjust epsilon) might be prohibitively
> difficult/slow? Or are there actually any situations where you'd want to
> compare two floating point numbers directly with each other?

So prohibitively difficult as not to be computable in some cases.

> Thanks for taking the time to answer my "newbie questions" :)

My pleasure.

Andrew.

Re: Floating point behavior

[Brian Dessent]

eool wrote:

> Basically, the reason given in the FAQ is that the result of cos(x)
> (intermediate, higher precision value) will be stored in a register and thus
> truncated. Then the result of cos(y) is directly compared with the truncated
> value located in a register, resulting in this strange (?) behavior. Now, my
> question is, why isn't gcc able to tell that it is about to compare two
> floating point numbers, one of which has been truncated, and then similarly

You can get this with -ffloat-store.  It's not enabled by default as it
decreases performance because it adds the redundant inefficiency of
having to transfer all the intermediate results of computations from
registers to memory only to just reload them back into registers again. 
Normally the optimizing compiler strives very hard to avoid this very
thing, so this is really a poor workaround to deal with the a quirk of
one particular architecture.

And in general floating point units do not work this way.  It's only the
x87 that has this quirk of performing calculations with extended
precision.  And even with this oddball implementation, the behavior is
still optional: the precision is configurable.  It's just that the
default is 80 bits on many operating systems because there are legacy
math routines that depend on the excess precision for numerical
stability.  But if you aren't dealing with such code you can simply set
the precision in the 387 control word to 64 bits and the problem
vanishes because there is no more 80->64 truncation.  Recent versions of
gcc even have an option (-mpc64) to do this for you at program startup.

An alternative is to avoid the 387 unit and its antiquated 21 year old
warts entirely.  The SSE unit has much more modern and sane design, and
it does calculations on doubles only in double precision.  You can use
-fpmath=sse to tell gcc to avoid the 387 unit and use the sse unit
instead.  Of course this means you have to target an architecture that
supports sse2 (for doubles) and that rules out a fair number of older
systems, notably the AMD line prior to the 64 bit models.  If you're
distributing software in binary form it may not be realistic to assume
sse2 support.  But for x86_64 targets, sse2 is required and is the
default anyway so you don't even have to really worry about it there.

> truncate the second value before the comparison? (In fact, since the return
> type of cos() is double, shouldn't the intermediate value be truncated to
> fit a double in any case?)

The IA-32 calling convention dictates that floating point return values
are passed in the %st(0) register, not on the stack.  So just calling a
function that returns double doesn't necessarily result in any
truncation.  It's only when a value is moved out of a register and
stored into memory that the truncation occurs.

> Second (this question might be even more naive than the previous one), if
> you are never, ever, supposed to compare floating point numbers using ==; if
> you are always supposed to check the difference of the numbers against an
> epsilon value; then how come the compiler can't automatically do this? Is it
> because analyzing the source code to track the accumulated floating point
> errors (to appropriately adjust epsilon) might be prohibitively
> difficult/slow? Or are there actually any situations where you'd want to
> compare two floating point numbers directly with each other?

C is not a language where the compiler has any room for interpretation
of what you specify.  People expect C compilers to generate the code
that's written without second guessing the intent of the programmer. 
Besides, the epsilon value has to be scaled appropriately to the
predicates of the comparison, so it's not like there's one value that
fits all situations.  And as already mentioned, there are existing
workarounds that mitigate or completely avoid the problem.

Again, the excess precision issue is a design quirk of a particular
implementation that, while at one point in time was very widespread, is
mercifully obsolete these days.  Users targeting sparc, power, arm,
mips, alpha (etc.) platforms would all be very unhappy to see their code
changed to work around something that doesn't affect them.

Brian

Re: Floating point behavior

[Vincent Lefevre]

On 2008-09-04 02:00:59 -0700, Brian Dessent wrote:
> eool wrote:
> > Basically, the reason given in the FAQ is that the result of
> > cos(x) (intermediate, higher precision value) will be stored in a
> > register and thus truncated. Then the result of cos(y) is directly
> > compared with the truncated value located in a register, resulting
> > in this strange (?) behavior.

There could be other reasons in practice, e.g. the cos() function
wasn't the same one. For instance, this can happen when one of them
is evaluated at compilation time because the argument value was
known.

> > Now, my question is, why isn't gcc
> > able to tell that it is about to compare two floating point
> > numbers, one of which has been truncated, and then similarly

Because it can't. You can have more complex expressions...

> You can get this with -ffloat-store.

But it has been told that it doesn't always have effects inside
expressions. You also need to store every intermediate value to
variables.

> And in general floating point units do not work this way. It's only
> the x87 that has this quirk of performing calculations with extended
> precision.

This is more a problem with operating systems (and API's) that enable
this extended precision. With some codes, there could still be a
problem with the extended exponent range, though.

>  And even with this oddball implementation, the behavior is
> still optional: the precision is configurable.  It's just that the
> default is 80 bits on many operating systems because there are legacy
> math routines that depend on the excess precision for numerical
> stability.

Only for math routines *specifically* written for the x87 configured
in extended precision, i.e. non-portable across the processors and
the operating systems, and not following the IEEE-754 standard (I
mean here that the code won't work on all IEEE-754 implementations
as extended precision, though allowed by the IEEE-754 standard, is
just optional). So, this is a completely stupid argument for having
extended precision by default.

-- 
Vincent Lefèvre <vincent@vinc17.org> - Web: <http://www.vinc17.org/>
100% accessible validated (X)HTML - Blog: <http://www.vinc17.org/blog/>
Work: CR INRIA - computer arithmetic / Arenaire project (LIP, ENS-Lyon)

Re: Floating point behavior

[Martin Krischik]

-----BEGIN PGP SIGNED MESSAGE-----
Hash: SHA1


Am 07.09.2008 um 12:00 schrieb Vincent Lefevre:

> On 2008-09-04 02:00:59 -0700, Brian Dessent wrote:
>> eool wrote:
>>> Basically, the reason given in the FAQ is that the result of
>>> cos(x) (intermediate, higher precision value) will be stored in a
>>> register and thus truncated. Then the result of cos(y) is directly
>>> compared with the truncated value located in a register, resulting
>>> in this strange (?) behavior.
>
> There could be other reasons in practice, e.g. the cos() function
> wasn't the same one. For instance, this can happen when one of them
> is evaluated at compilation time because the argument value was
> known.
>

Indeed - all trigonometric functions are only aproximations. I only  
recently had a good look into it.

http://uiq3.sourceforge.net/wiki/index.php/Office/FX-602P/Diff#Arithmetic_Libraries

Martin
- --
Martin Krischik
krischik@users.sourceforge.net

-----BEGIN PGP SIGNATURE-----
Version: GnuPG v1.4.8 (Darwin)

iD8DBQFIw6xqijwKaHyem9cRAtLoAJ4g+aPY/QD5Wz+tGRQ+HSI9W9kOEwCbBC7u
o6rTCrX9wMuv3nm27U4r2P8=
=1opA
-----END PGP SIGNATURE-----

Re: Floating point behavior

[Vincent Lefevre]

On 2008-09-07 12:26:50 +0200, Martin Krischik wrote:
> Indeed - all trigonometric functions are only aproximations. I only  
> recently had a good look into it.
>
> http://uiq3.sourceforge.net/wiki/index.php/Office/FX-602P/Diff#Arithmetic_Libraries

You can also look at:

  http://www.vinc17.org/research/testlibm/

-- 
Vincent Lefèvre <vincent@vinc17.org> - Web: <http://www.vinc17.org/>
100% accessible validated (X)HTML - Blog: <http://www.vinc17.org/blog/>
Work: CR INRIA - computer arithmetic / Arenaire project (LIP, ENS-Lyon)