Peculiar floating point numbers in GCC

[n.torrey.pines]

I understand that with floats, x*y == y*x, for example, might not
hold. But what if it's the exact same operation on both sides?

I tested this with GCC 3.4.4 (Cygwin), and it prints 0 (compiled with -
g flag)

#include <cmath>
#include <iostream>

int main() {
double x = 3.0;

double f = std::cos(x);
std::cout << (f == std::cos(x)) << std::endl;
}

Using -O3 instead optimizes something, and now both sides are equal.

VC++8.0 prints 1 in both debug and release.

What does the standard have to say?

 

[Jim Langston]

I understand that with floats, x*y == y*x, for example, might not
hold. But what if it's the exact same operation on both sides?

I tested this with GCC 3.4.4 (Cygwin), and it prints 0 (compiled with -
g flag)

#include <cmath>
#include <iostream>

int main() {
double x = 3.0;

double f = std::cos(x);
std::cout << (f == std::cos(x)) << std::endl;
}

Using -O3 instead optimizes something, and now both sides are equal.

VC++8.0 prints 1 in both debug and release.

What does the standard have to say?

It's not really the standard that's the issue I dont' think, it's just the
way floating point math works.

In your particular case, those statements are close together, initializing f
and the comparison, so the compiler may be optimizing and comparing the same
thing. It all depends on how the compiler optimizes. Even something like:

double f = std::cos(x);
double g = std::cos(x);
std::cout << ( f == g ) << std::endl;

may output 1 or 0, depending on compiler optimization.

You just cant count on floating point equality, it may work sometimes, not
others.

 

[Pete Becker]

It's not really the standard that's the issue I dont' think, it's just the
way floating point math works.

In your particular case, those statements are close together, initializing f
and the comparison, so the compiler may be optimizing and comparing the same
thing. It all depends on how the compiler optimizes. Even something like:

double f = std::cos(x);
double g = std::cos(x);
std::cout << ( f == g ) << std::endl;

may output 1 or 0, depending on compiler optimization.

It had better output 1, regardless of compiler optimizations.

You just cant count on floating point equality, it may work sometimes, not
others.

The not-so-subtle issue here is that the compiler is permitted to do
floating-point arithmetic at higher precision than the type requires. On
the x86 this means that floating-point math is done with 80-bit values,
while float is 32 bits and double is 64 bits. (The reason for allowing
this is speed: x86 80-bit floating-point math is much faster than 64-bit
math) When you store into a double, the stored value has to have the
right precision, though. So storing the result of cos in a double can
force truncation of a wider-than-double value. When comparing the stored
value to the original one, the compiler widens the original out to 80
bits, and the result is different from the original value. That's what's
happening in the original example.

In your example, though, both results are stored, so when widened they
should produce the same result. Some compilers don't do this, though,
unless you tell them that they have to obey the rules (speed again).

--

-- Pete
Roundhouse Consulting, Ltd. (www.versatilecoding.com)
Author of "The Standard C++ Library Extensions: a Tutorial and
Reference." (www.petebecker.com/tr1book)

 

[Jim Langston]

It had better output 1, regardless of compiler optimizations.


The not-so-subtle issue here is that the compiler is permitted to do
floating-point arithmetic at higher precision than the type requires. On
the x86 this means that floating-point math is done with 80-bit values,
while float is 32 bits and double is 64 bits. (The reason for allowing
this is speed: x86 80-bit floating-point math is much faster than 64-bit
math) When you store into a double, the stored value has to have the right
precision, though. So storing the result of cos in a double can force
truncation of a wider-than-double value. When comparing the stored value
to the original one, the compiler widens the original out to 80 bits, and
the result is different from the original value. That's what's happening
in the original example.

In your example, though, both results are stored, so when widened they
should produce the same result. Some compilers don't do this, though,
unless you tell them that they have to obey the rules (speed again).

FAQ 29.18 disagrees with you.

http://www.parashift.com/c++-faq-lite/newbie.html#faq-29.18

 

[James Kanze]

Jim Langston wrote:

Don't get paranoid. <g> There's a specific reason for this descrepancy.

But it is partially the standard at fault, since the standard
explicitly allows extended precision for intermediate results
(including, I think function return values, but I'm not sure
there).

It had better output 1, regardless of compiler optimizations.

I don't know about this particular case, but I've had similar
cases fail with g++. By default, g++ is even looser than the
standard allows. One could put forward the usual argument about
it not mattering how fast you get the wrong results, but in at
least some cases, the results that it gets by violating the
standard are also more precise. (The problem is, of course,
that precise or not, they're not very predictable.)

The not-so-subtle issue here is that the compiler is permitted to do
floating-point arithmetic at higher precision than the type requires. On
the x86 this means that floating-point math is done with 80-bit values,
while float is 32 bits and double is 64 bits. (The reason for allowing
this is speed: x86 80-bit floating-point math is much faster than 64-bit
math) When you store into a double, the stored value has to have the
right precision, though. So storing the result of cos in a double can
force truncation of a wider-than-double value. When comparing the stored
value to the original one, the compiler widens the original out to 80
bits, and the result is different from the original value. That's what's
happening in the original example.

Question: if I have a function declared "double f()", is the
compiler also allowed to maintain extended precision, even
though I'm using the results of the expression in the return
statement to initialize a double? (G++ definitly does, at least
in some cases.)

In your example, though, both results are stored, so when widened they
should produce the same result. Some compilers don't do this, though,
unless you tell them that they have to obey the rules (speed again).

OK. We're on the same wave length. I can confirm that g++ is
one of those compilers.

--

 

[Walter Bright]

Question: if I have a function declared "double f()", is the
compiler also allowed to maintain extended precision, even
though I'm using the results of the expression in the return
statement to initialize a double? (G++ definitly does, at least
in some cases.)

I don't quite understand your question. I would rephrase it as "if a
function is declared to return a double, must it round the return value
to double precision, or can it return a value that actually is of higher
precision?"

For the Digital Mars C/C++, the answer is "no". Even though the return
value is typed as a double, it is actually returned as a full 80 bit
value in the ST0 register. In other words, it is treated as an
intermediate value that can be of higher precision than the type.

D programming is even more liberal with floating point - the compiler is
allowed to do all compile time constant folding, even copy propagation,
at the highest precision available.

The idea is that one should structure algorithms so that they work with
the minimum precision specified by the type, but don't break if a higher
precision is used.

 

[Walter Bright]

For the Digital Mars C/C++, the answer is "no".

Er, I meant that it returns a value that is actually of higher precision.

 

[Bo Persson]

Jim Langston wrote:
:::
::: In your example, though, both results are stored, so when widened
::: they should produce the same result. Some compilers don't do this,
::: though, unless you tell them that they have to obey the rules
::: (speed again).
::
:: FAQ 29.18 disagrees with you.
::
:: http://www.parashift.com/c++-faq-lite/newbie.html#faq-29.18

The FAQ tells you what actually happens when you don't tell the compiler to
obey the language rules.

Like Pete says, most compilers have an options to select fast or strictly
correct code. The default is of course fast, but slightly incorrect.

If the source code stores into variables x and y, the machine code must do
that as well. Unless compiler switches (or lack thereof) relax the
requirements.


The "as-if" rule of code transformations doesn't apply here, as we can
actually see the difference.


Bo Persson

 

[James Kanze]

Jim Langston wrote:

:: "Pete Becker" <[email protected]> wrote in message

::: In your example, though, both results are stored, so when widened
::: they should produce the same result. Some compilers don't do this,
::: though, unless you tell them that they have to obey the rules
::: (speed again).

:: FAQ 29.18 disagrees with you.

The FAQ tells you what actually happens when you don't tell the compiler to
obey the language rules.

Like Pete says, most compilers have an options to select fast or strictly
correct code. The default is of course fast, but slightly incorrect.

For a one definition of "correct": for a naive definition,
1.0/3.0 will never give correct results, regardless of
conformance.

And I don't understand "slightly"---I always thought that
"correct" was binary: the results are either correct, or they
are not.

If the source code stores into variables x and y, the machine code must do
that as well.

Just a nit: the results must be exactly the same as if the
machine code had stored into the variables. The memory write is
not necessarily necessary, but any rounding that would have
occured is.

Unless compiler switches (or lack thereof) relax the
requirements.

The "as-if" rule of code transformations doesn't apply here, as we can
actually see the difference.

The "as-if" rule always applies, since it says that only the
results count. In this case, for example: I believe that the
Intel FPU has an instruction to force rounding to double
precision, without actually storing the value, and an
implementation could use this.

--

 

[James Kanze]

I don't quite understand your question. I would rephrase it as "if a
function is declared to return a double, must it round the return value
to double precision, or can it return a value that actually is of higher
precision?"

Yes. Your English is better than mine.

For the Digital Mars C/C++, the answer is "no". Even though the return
value is typed as a double, it is actually returned as a full 80 bit
value in the ST0 register. In other words, it is treated as an
intermediate value that can be of higher precision than the type.

My question concerned what the standard requires, not what
different compilers do. I know what the compilers I use do (and
that it differs according to the platforms); my question perhaps
belongs more in comp.std.c++, but I wanted to know whether they
were conform, or whether it was "yet another bug" in the
interest of getting the wrong results faster.

D programming is even more liberal with floating point - the compiler is
allowed to do all compile time constant folding, even copy propagation,
at the highest precision available.

Which has both advantages and disadvantages. As I pointed out
elsewhere, the "non-conformant" behavior of g++ can sometimes
result in an answer that is more precise that conformant
behavior would have been. (Extended precision was invented for
a reason.)

The idea is that one should structure algorithms so that they work with
the minimum precision specified by the type, but don't break if a higher
precision is used.

The problem is that if you want to analyse your limits in
detail, you don't really know what precision is being used.
There are two schools of thought here, and I don't know enough
numerics to have an opinion as to which is right. (The only
calculations in my present application involve monetary amounts,
so still another set of rules applies.)

 

[Pete Becker]

FAQ 29.18 disagrees with you.

http://www.parashift.com/c++-faq-lite/newbie.html#faq-29.18

No, it says the same thing, as far as it goes. It doesn't mention the
requirement that the compiler respect stored values, and, as I said,
many don't do that unless you explicitly ask for it.

--

-- Pete
Roundhouse Consulting, Ltd. (www.versatilecoding.com)
Author of "The Standard C++ Library Extensions: a Tutorial and
Reference." (www.petebecker.com/tr1book)

 

[Walter Bright]

The problem is that if you want to analyse your limits in
detail, you don't really know what precision is being used.
There are two schools of thought here, and I don't know enough
numerics to have an opinion as to which is right. (The only
calculations in my present application involve monetary amounts,
so still another set of rules applies.)

My experience in the matter is that the only algorithms that failed when
used with greater precision were:

1) test suites
2) wrong for other reasons, too

Put another way, if I was implementing a square root algorithm, why
would I ever *need* less accuracy?

Ok, ok, I can think of one - I'm writing an emulator that needs to be
dumbed down to match its dumb target.

 

[James Kanze]

My experience in the matter is that the only algorithms that failed when
used with greater precision were:

1) test suites
2) wrong for other reasons, too

Put another way, if I was implementing a square root algorithm, why
would I ever *need* less accuracy?

I'm sure that that is true for many algorithms. But the
question isn't more or less accuracy, it is known accuracy, and
I can imagine certain algorithms becoming instable if the
accuracy varies.

As I said, I don't know the domain enough to have an opinion. I
do know that in my discussions with experts in the domain,
opinions varied, and many seem to consider the extended accuracy
in the Intel processors a mis-design (or maybe I've
misunderstood them). Curiously, all seem to consider the base
16 used in IBM mainframe floating point bad thing, and the most
common reason they give is the varying accuracy.

 

[n.torrey.pines]

http://www.parashift.com/c++-faq-lite/newbie.html#faq-29.18

Thanks, that was very useful. I like the tone of the FAQ though: "if
you don't think this is surprising, you are stupid, but if you do,
that's because you don't know very much - keep reading"

 

[Walter Bright]

In this case, for example: I believe that the
Intel FPU has an instruction to force rounding to double
precision, without actually storing the value, and an
implementation could use this.

That rounds the mantissa, but not the exponent. Thus, it doesn't do the
job. The only way to get the 87 to create a true double value is to
store it into memory, a slow operation.

Early Java's insistence on doubles for intermediate values caused a lot
of problems because of this.

 

[Bo Persson]

James Kanze wrote:
:: Jim Langston wrote:
:
:
::::: In your example, though, both results are stored, so when widened
::::: they should produce the same result. Some compilers don't do this,
::::: though, unless you tell them that they have to obey the rules
::::: (speed again).
:
:::: FAQ 29.18 disagrees with you.
:
:::: http://www.parashift.com/c++-faq-lite/newbie.html#faq-29.18
:
:: The FAQ tells you what actually happens when you don't tell the
:: compiler to obey the language rules.
:
:: Like Pete says, most compilers have an options to select fast or
:: strictly correct code. The default is of course fast, but slightly
:: incorrect.
:
: For a one definition of "correct": for a naive definition,
: 1.0/3.0 will never give correct results, regardless of
: conformance.

It has an expected result, given the hardware. Like you mentioned elsewhere,
we have a problem with the Intel hardware, that a temporary result (in a
register) is produced faster *and* has higher precision than a stored
result.

That makes it very tempting to use the temporary.

: And I don't understand "slightly"---I always thought that
: "correct" was binary: the results are either correct, or they
: are not.

That was supposed to be some kind of humor. Didn't work out too well, did
it?

:
:: If the source code stores into variables x and y, the machine code
:: must do that as well.
:
: Just a nit: the results must be exactly the same as if the
: machine code had stored into the variables. The memory write is
: not necessarily necessary, but any rounding that would have
: occured is.
:
:: Unless compiler switches (or lack thereof) relax the
:: requirements.
:
:: The "as-if" rule of code transformations doesn't apply here, as we
:: can actually see the difference.
:
: The "as-if" rule always applies, since it says that only the
: results count. In this case, for example: I believe that the
: Intel FPU has an instruction to force rounding to double
: precision, without actually storing the value, and an
: implementation could use this.

You must do a store to get the rounding. However, the target of the store
instruction can be another FP register, so in effect Yes.

In the example from the FAQ, the machine code compares a stored value to a
temporary in a register. That doesn't follow the "as-if" rule (as we can
easily see). Possibly because the code wasn't compiled with a "strict but
slow" compiler option.


Bo Persson

 

[James Kanze]

That rounds the mantissa, but not the exponent. Thus, it doesn't do the
job.

I'm not sure what you mean by not rounding the exponent. The
exponent is an integral value, so no rounding is involved.

If it doesn't truncate the exponent, that's not a problem, since
anytime it would have to truncate the exponent, you have
undefined behavior. According to the standard, anyway; from a
quality of implementation point of view, I don't know. Most
implementations (those using IEEE, anyway) *do* define behavior
in this case.

The only way to get the 87 to create a true double value is to
store it into memory, a slow operation.

Even rounding in a register takes measurable time. (But not as
much as storing it, of course.)

Early Java's insistence on doubles for intermediate values caused a lot
of problems because of this.

I know. Their goal was laudable: that code should always give
the same results, everywhere. The problem is that in order to
do so, floating point becomes real slow on some platforms,
including the Intel. The other alternative would have been
requiring extended precision in the intermediate values, but
then, they couldn't use native floating point on a Sparc, which
doesn't support extended precision. And for some reason, Sun
seems to consider Sparcs an important platform. So they
compromized. (The still require IEEE, which makes Java
significantly slower on an IBM mainframe. Also an important
platform in certain milieu. Obviously, two measures apply.)

 

[James Kanze]

James Kanze wrote:

: On Apr 7, 9:08 am, "Bo Persson" <[email protected]> wrote::: Jim Langston wrote:

:::: "Pete Becker" <[email protected]> wrote in message

::::: In your example, though, both results are stored, so when widened
::::: they should produce the same result. Some compilers don't do this,
::::: though, unless you tell them that they have to obey the rules
::::: (speed again).

:::: FAQ 29.18 disagrees with you.

:: The FAQ tells you what actually happens when you don't tell the
:: compiler to obey the language rules.

:: Like Pete says, most compilers have an options to select fast or
:: strictly correct code. The default is of course fast, but slightly
:: incorrect.

: For a one definition of "correct": for a naive definition,
: 1.0/3.0 will never give correct results, regardless of
: conformance.

It has an expected result, given the hardware. Like you mentioned elsewhere,
we have a problem with the Intel hardware, that a temporary result (in a
register) is produced faster *and* has higher precision than a stored
result.

That's why I raised the question of what is meant by correct.
In context, if I remember it right, Pete's statement was clear:
correct meant conform to the standards. The simple quote above
removed the context, however, and makes the statement somewhat
ambiguous. Correct can mean many things, especially where
floating point values are involved.

That makes it very tempting to use the temporary.

Quite. Especially because in some (most?, all?) cases, the
resulting code will be "correct" in the sense that it meets its
requirement specifications.

I think that this is the basis of Walter's argument. He hasn't
convinced me that it applies in all cases, but even before this
thread, I was convinced that it applied often enough that a
compiler should offer the alternative. And while I'd normally
oppose non-standard as a default, floating point is so subtle
that either way, someone who doesn't understand the issues will
get it wrong, so the argument of "naïve users getting the
correct behavior" really doesn't apply.

: And I don't understand "slightly"---I always thought that
: "correct" was binary: the results are either correct, or they
: are not.

That was supposed to be some kind of humor. Didn't work out too well, did
it?

No humor. For any given definition of "correct", code is either
correct, or it is not correct.

:: If the source code stores into variables x and y, the machine code
:: must do that as well.

: Just a nit: the results must be exactly the same as if the
: machine code had stored into the variables. The memory write is
: not necessarily necessary, but any rounding that would have
: occured is.

:: Unless compiler switches (or lack thereof) relax the
:: requirements.

:: The "as-if" rule of code transformations doesn't apply here, as we
:: can actually see the difference.

: The "as-if" rule always applies, since it says that only the
: results count. In this case, for example: I believe that the
: Intel FPU has an instruction to force rounding to double
: precision, without actually storing the value, and an
: implementation could use this.

You must do a store to get the rounding. However, the target of the store
instruction can be another FP register, so in effect Yes.

I didn't remember it that way, but it's been years since I
looked at the specification. (It was the specification for the
original 8087, so that gives you some idea of just how long ago
it was.)

In the example from the FAQ, the machine code compares a stored value to a
temporary in a register. That doesn't follow the "as-if" rule (as we can
easily see). Possibly because the code wasn't compiled with a "strict but
slow" compiler option.

Quite possibly. I know that real programs do occasionally have
this problem. I seem to recall that g++ even has it in certain
cases where the strict option is given, but I could be wrong.
This problem is probably the strongest argument against extended
precision. The fact that the value may change depending on
whether an intermediate value was spilled to memory or
maintained in a register certainly cannot help reasoning about
program correctness.

 

[Walter Bright]

I'm not sure what you mean by not rounding the exponent. The
exponent is an integral value, so no rounding is involved.

I meant that the exponent for 80 bit doubles has a greater range, and
values in that greater range do not cause the whole to be set to
infinity when set to double precision significand.