floating point arithmetic

[Szabolcs Nagy]

hello, i get unexpected floating point behaviour

the following code prints "not eq" on various compilers with various
settings on my 32bit x86 machine:

#include <stdlib.h>
#include <stdio.h>

int main() {
if ((double)atoi("1")/atoi("3") == (double)atoi("1")/atoi
("3"))
puts("eq");
else
puts("not eq");
return 0;
}

with gcc it seems the first division is performed with x86 80bit float
registers but the second division is done with a different instruction
which i don't know exactly

i know floating point equality check is rarely adequate and there are
rounding errors etc.

but does the standard allow such behaviour in this case?
ie. is it a compiler bug or 'undefined behaviour'?

 

[Tim Prince]

hello, i get unexpected floating point behaviour

the following code prints "not eq" on various compilers with various
settings on my 32bit x86 machine:

#include <stdlib.h>
#include <stdio.h>

int main() {
if ((double)atoi("1")/atoi("3") == (double)atoi("1")/atoi
("3"))
puts("eq");
else
puts("not eq");
return 0;
}

with gcc it seems the first division is performed with x86 80bit float
registers but the second division is done with a different instruction
which i don't know exactly

This depends on which options you set, and which target you compile for.
Please read compiler manuals, if questions, ask on the associated help forum.
Also, some constant propagation calculations are done at compile time,
others at run time.
In a case like this, if you chose, you could make a stronger effort to
have them all done the same.

i know floating point equality check is rarely adequate and there are
rounding errors etc.

but does the standard allow such behaviour in this case?
ie. is it a compiler bug or 'undefined behaviour'?

In my Harbison & Steele it says, in effect, that agreement within 1 ULP is
the best to be expected.

 

[Keith Thompson]

hello, i get unexpected floating point behaviour

the following code prints "not eq" on various compilers with various
settings on my 32bit x86 machine:

#include <stdlib.h>
#include <stdio.h>

int main() {
if ((double)atoi("1")/atoi("3") == (double)atoi("1")/atoi
("3"))
puts("eq");
else
puts("not eq");
return 0;
} [...]
but does the standard allow such behaviour in this case?
ie. is it a compiler bug or 'undefined behaviour'?

It's not undefined behavior, it's merely unspecified. The program
will print either "eq" or "not eq".

("int main()" is better written as "int main(void)".)

 

[Szabolcs Nagy]

It's not undefined behavior, it's merely unspecified. The program
will print either "eq" or "not eq".

i'm not a language lawyer but there are certain requirements by the
standard and it seemed suspicious to me that the division (and other
arithmetic operators) should be well defined as a requirement either
in c or in ieee 754 (which is required by c99 afaik)

where well defined mean that the same operation on the same numbers of
the same type always gives the same result on a given architecture/
compiler/..

anyway thanks for your reply

 

[Ben Bacarisse]

i'm not a language lawyer but there are certain requirements by the
standard and it seemed suspicious to me that the division (and other
arithmetic operators) should be well defined as a requirement either
in c or in ieee 754 (which is required by c99 afaik)

C99 does not require ieee 754 floating point (BTW I will call it IFP
because the alternative is both too long and the new IEC number is
less memorable). However where the implementation is based on IFP, as
here, the C standard essentially defers to it. However, the only
requirement from the FP standard is on the minimum accuracy of the
result. C goes on to permit (but does require) temporary calculations
using more accuracy if available.

It is likely that which ever calculation is done second is done in 80
bits so that the registers involved have the maximum precision
available. If the earlier result is loaded from memory, it will be a
mere 64 bits. It is possible the compiler knew the result was to be
stored and hence it may have chosen a shorter divide operation. All
that is speculative -- I know very little about Intel FP -- but the
possibility that the same calculation does not compare equal to itself
is specifically permitted by C, provided both values are as accurate
as the IFP standard mandates.

The compare not permitted to be not equal because one of the two is
not accurate enough; but one may have too many bits to compare equal
to the other.

where well defined mean that the same operation on the same numbers of
the same type always gives the same result on a given architecture/
compiler/..

In short, the resul must always be good enough but not always the same
even on a single machine.

 

[jameskuyper]

i'm not a language lawyer but there are certain requirements by the
standard and it seemed suspicious to me that the division (and other
arithmetic operators) should be well defined ...

You're right to be suspicious of that; it's not actually the case.

...as a requirement either
in c or in ieee 754 (which is required by c99 afaik)

I believe that IEEE 754 is not mentioned in C90; it is not mandatory
in C99. If an implementation pre#defines __STDC_IEC_559_ it is
supposed to conform to IEC/IEEE 60559, which is equivalent to IEEE
754, but otherwise it's free to do whatever it wants.

I don't have a copy of IEEE 754; my copy of IEEE 854 (which is a radix-
independent standard that is otherwise similar to IEEE 754) has gone
missing. However, IIRC, even if __STDC_IEE_559__ is pre#defined,
differences of 1ulp are allowed.

where well defined mean that the same operation on the same numbers of
the same type always gives the same result on a given architecture/
compiler/..

'undefined behavior' is a piece of C jargon that is very precisely
defined by the C standard to have a meaning that is not equivalent to
"not well-defined". Undefined behavior means, in the context of the C
standard, behavior on which the C standard imposes no requirements.
It's quite possible (and in many cases likely) that some other source
(the POSIX standard, a platform-specific ABI standard, the User's
guide for a given implementation of C) does in fact provide a
definition; but it's still "undefined behavior" as far as the C
standard is concerned. The C standard might define the behavior
incompletely, in which case the behavior is "unspecified", not
"undefined".

In this case, the C standard does require that the result of the
division must be, "the quotient from the division of the first operand
by the second" (6.5.5p5), so the behavior is not undefined. However,
the C standard does not specify the accuracy with which that quotient
is calculated, so the result is unspecified..

 

[bartc]

hello, i get unexpected floating point behaviour

the following code prints "not eq" on various compilers with various
settings on my 32bit x86 machine:

#include <stdlib.h>
#include <stdio.h>

int main() {
if ((double)atoi("1")/atoi("3") == (double)atoi("1")/atoi
("3"))
puts("eq");
else
puts("not eq");
return 0;
}

with gcc it seems the first division is performed with x86 80bit float
registers but the second division is done with a different instruction
which i don't know exactly

i know floating point equality check is rarely adequate and there are
rounding errors etc.

but does the standard allow such behaviour in this case?
ie. is it a compiler bug or 'undefined behaviour'?

It does seem to be a 'bug' due to the mixed precision of registers and
memory, as I think Ben mentioned;

I tried an extra (double) cast around each of the expressions, and would
expected this to have reduced both expressions down to 64-bits, but perhaps
it doesn't bother when the expression is already in an 80-bit fp register. I
don't think it's that easy, or efficient, to do this on x86 FPU. Or it
decided the result was already 'double' even though one side might really
have been 'long double' at that point.

FWIW, lcc-win32 and DMC compilers for x86 both give an "eq" result on your
code.

 

[Keith Thompson]

C99 does not require ieee 754 floating point (BTW I will call it IFP
because the alternative is both too long and the new IEC number is
less memorable). However where the implementation is based on IFP, as
here, the C standard essentially defers to it. However, the only
requirement from the FP standard is on the minimum accuracy of the
result. C goes on to permit (but does require) temporary calculations
using more accuracy if available.

[...]

Typo: you meant "(but does not require)".

 

[Szabolcs Nagy]

thanks for the replies, both were informational

In this case, the C standard does require that the result of the
division must be, "the quotient from the division of the first operand
by the second" (6.5.5p5), so the behavior is not undefined. However,
the C standard does not specify the accuracy with which that quotient
is calculated, so the result is unspecified..

you (and the standard) mention accuracy, however accuracy is not the
issue here

it does not matter if (double)1/3 evaluates to 0.3 or 0.4, what
matters whether the expression is idempotent or not (always give the
same answer)

from a theoretical point of view idempotency is very important, if it
is not present then calculations are not reproducible, the program is
not observable (adding a print statement might change the result) etc.

from implementation or practical point of view requiring idempotency
might be too much, at least it seems that's what the standardisation
committee thought

 

[Ben Bacarisse]

However, the only
requirement from the FP standard is on the minimum accuracy of the
result. C goes on to permit (but does require) temporary calculations
using more accuracy if available.

[...]

Typo: you meant "(but does not require)".

*sigh* Yes. Thank you.

 

[James Kuyper]

thanks for the replies, both were informational


you (and the standard) mention accuracy, however accuracy is not the
issue here

it does not matter if (double)1/3 evaluates to 0.3 or 0.4, what
matters whether the expression is idempotent or not (always give the
same answer)

The phrase you're looking for is probably "repeatable" or "consistent",
not "idempotent". When applied to a binary operation like division,
"idempotent" means that the result of applying the operator to equal
operands always produces a result equal to either operand. For instance,
maximum(a,b) is idempotent because maximum(a,a)==a. Division is most
certainly not idempotent; the only value of x for which x/x == x is 1.

The standard imposes no requirements on the accuracy; in particular, it
imposes no requirement that quotients involving identical operands be
evaluated with the same accuracy every time they are evaluated.

from a theoretical point of view idempotency is very important, if it
is not present then calculations are not reproducible, the program is
not observable (adding a print statement might change the result) etc.

Consistency is a useful property to have in general, and I can't think
of any reason why an implementation would want to have this kind of
inconsistency. However, floating point operations are inherently inexact
in general, and any code which is written with the expectation of that
mathematically equivalent expressions should be numerically exactly
equal is broken at a deep conceptual level.

 

[Szabolcs Nagy]

The phrase you're looking for is probably "repeatable" or "consistent",
not "idempotent". When applied to a binary operation like division,
"idempotent" means that the result of applying the operator to equal
operands always produces a result equal to either operand. For instance,
maximum(a,b) is idempotent because maximum(a,a)==a. Division is most
certainly not idempotent; the only value of x for which x/x == x is 1.

well you are right, however in computer science idempotence is used
with a slightly different meaning:
http://en.wikipedia.org/wiki/Idempotence#In_computing
now i see "referentially transparent" or "pure function" would be a
better term for it

The standard imposes no requirements on the accuracy; in particular, it
imposes no requirement that quotients involving identical operands be
evaluated with the same accuracy every time they are evaluated.

i meant that it mentions accuracy like "accuracy of operations are
implementation specific" but not that "the operations might not be
consistent"

in general, and any code which is written with the expectation of that
mathematically equivalent expressions should be numerically exactly
equal is broken at a deep conceptual level.

ok i can accept that, but there are cases when it would make life a
bit easier:
assume a numerically unstable simulation, something happens at the
1000th step that should not happen, the programmer binds some
subexpression of an internal calculation to a variable to debug the
code, which changes the outcome of the program, so the problem is
gone, the program is not debuggable.

 

[Flash Gordon]

James Kuyper wrote:


ok i can accept that, but there are cases when it would make life a
bit easier:
assume a numerically unstable simulation, something happens at the
1000th step that should not happen, the programmer binds some
subexpression of an internal calculation to a variable to debug the
code, which changes the outcome of the program, so the problem is
gone, the program is not debuggable.

It is one instance where I would use a debugger. A decent debugger will
not interfere with the values (I've come across hardware
in-circuit-emulators and simulators where things did miss-behave, but
that is the exception not the rule).

In any case, there are good reasons compilers do not always do the same
thing. The processor might be able to do more things in parallel if it
selects a different instruction, or it may save taking things out of
registers in to RAM, or... the possible reasons are endless. This is
also why C99 introduced pragmas which allow you to tell the compiler to
*not* do some of these things, they are so valuable they are needed but
sometimes one needs the maths to be more well behaved.

 

[Phil Carmody]

The phrase you're looking for is probably "repeatable" or
"consistent", not "idempotent". When applied to a binary operation
like division, "idempotent" means that the result of applying the
operator to equal operands always produces a result equal to either
operand.

Hmm, not really. It doesn't directly apply to a binary operator.
However, there are a couple of ways you can get to a unary operator
from a binary one, and then it can be applied. One way is duplication
of the single operand, which is what you describe, but the other way
is to curry the binary operator into a unary one by using one of the
operands.

The confusion between the two probably because square matrices can
act upon themselves by multiplication (thus represent both the operator
and the single operand to that operator). Viewing a matrix as an
operator is effectively currying the multiplication.

Phil

 

[James Kuyper]

James Kuyper wrote: ....

i meant that it mentions accuracy like "accuracy of operations are
implementation specific" but not that "the operations might not be
consistent"

If the standard imposed requirements on either the accuracy or the
consistency, it could in fact have chosen to separate the concepts. The
only thing it says is that the accuracy is implementation-defined, and
it goes so far as to say that the implementation-provided definition of
the accuracy can be "the accuracy is unknown". Given what the standard
fails to say about this issue, it would be entirely permissible for an
implementation to define that the accuracy is different depending upon
whether the operation is executed on an odd or even clock cycle. That's
just one of the most plausible of a great many inconvenient possibilities.

ok i can accept that, but there are cases when it would make life a
bit easier:
assume a numerically unstable simulation, something happens at the
1000th step that should not happen, the programmer binds some
subexpression of an internal calculation to a variable to debug the
code, which changes the outcome of the program, so the problem is
gone, the program is not debuggable.

Yes, that's pretty much what happens. You often have to turn off
optimizations to make any sense of the behavior of a program. I've used
debuggers that allow you to monitor the behavior of optimized code - but
it doesn't make any sense: executing a statement that is supposed to
update the value of a variable does not in fact do so, because the
optimizer has moved the update of the variable's value to another point
in the code. Turn off optimizations, and the bug you're looking for
often disappears. This is a problem with all kinds of code, but it's
more of a problem with floating point code, because the inherent
inaccuracy of floating point operations gives implementations more
freedom to implement optimizations that actually change the result. This
is the way real compilers work, and the standard was written to
accommodate that fact.

 

[Ben Bacarisse]

It does seem to be a 'bug' due to the mixed precision of registers and
memory, as I think Ben mentioned;

I see 'bug' in scare quotes but just to be clear, I don't see this as
a bug in gcc.

I tried an extra (double) cast around each of the expressions, and
would expected this to have reduced both expressions down to 64-bits,
but perhaps it doesn't bother when the expression is already in an
80-bit fp register. I don't think it's that easy, or efficient, to do
this on x86 FPU. Or it decided the result was already 'double' even
though one side might really have been 'long double' at that point.

FWIW, lcc-win32 and DMC compilers for x86 both give an "eq" result on
your code.

I am not sure that says very much. For example:

gcc 4.3.3 -O0 gives not eq
gcc 4.3.3 -O0 -ffloat-store gives eq
gcc 4.3.3 -O[123] give eq
icc 10.1 -O[0123] give eq

 

[Ben Bacarisse]

Puzzling.

#include <stdio.h>
#include <stdlib.h>
int main(void) {
double a, b;
if ((a=(double)atoi("1")/atoi("3")) ==
(b=(double)atoi("1")/atoi("3")))
puts("equal");
else
puts("not equal");
if (a == b)
puts("equal");
else
puts("not equal");
return 0;
}

Prints
not equal
equal

gcc 3.1

gcc 4.3.3

-O0 equal, equal
-O1 not equal, equal
-O1 -ffloat-store equal, equal

icc 10.0

-O0 equal, equal
-O1 not equal, no equal

You can get pretty much any result you like (although "equal" and "not
equal" will be less likely I think) but all are valid results in my
opinion.

 

[Tim Prince]

gcc 4.3.3

-O0 equal, equal
-O1 not equal, equal
-O1 -ffloat-store equal, equal

icc 10.0

-O0 equal, equal
-O1 not equal, no equal

You can get pretty much any result you like (although "equal" and "not
equal" will be less likely I think) but all are valid results in my
opinion.

gcc 4.5 win32
-O -march=pentium-m -mfpmath-sse equal,equal

icl 11.1 win32
-prec-div equal,equal

SSE2 platforms have been the norm for 8 years, and most compilers have
caught up.

 

[Mark Dickinson]

<snip>









gcc 4.3.3

-O0               equal,     equal
-O1               not equal, equal
-O1 -ffloat-store equal,     equal

icc 10.0

-O0               equal,     equal
-O1               not equal, no equal

You can get pretty much any result you like (although "equal" and "not
equal" will be less likely I think) but all are valid results in my
opinion.

Can someone explain why this behaviour doesn't contravene
C99 5.1.2.3, paragraph 12? "... In particular, casts and
assignments are required to perform their specified
conversion"? I realize that 5.1.2.3p12 is an example, and
hence non-normative, but it does still seem to show that this
behaviour is counter to the *intent* of the standard.

Maybe I'm misreading this paragraph, but to me it reads as
saying that in Joe's example, a and b should have been
rounded to double from whatever extended precision they
were computed in, *before* either of the equalities is applied.

 

[Ben Bacarisse]

gcc 4.5 win32
-O -march=pentium-m -mfpmath-sse equal,equal

icl 11.1 win32
-prec-div equal,equal

SSE2 platforms have been the norm for 8 years, and most compilers have
caught up.

The same applies to 4.3.3, of course, but I did not want to post a set
of unsurprising results!