double to int conversion yields strange results

[=?ISO-8859-1?Q?Bj=F8rn_Augestad?=]

Michael Mair wrote:




How does this differ from "n = 1.0 / d"?

1.0/d may have been computed with higher precision and the result
may be there in higher precision, too. By telling the compiler
explicitly that I want the result converted to double (before
it is converted to int), I would expect the excess precision
to disappear.
[snip]

In my understanding, inserting the cast should lead to the same
result as storing into an intermediate variable of the type we cast
to and using this variable for the next operation.

Maybe I understand something wrong but with the cast I would
expect the original example to work as expected by the OP on a
conforming implementation.

So would I, but a cast didn't help. I tried all kinds of casts and
tricks prior to posting the OP, but could not find a clean solution to
the problem. The code in the OP was just a simplified example of the
actual code, which all in all is close to 10KLOC with financial formulas.

Bjørn

 

[Christian Kandeler]

1.0/d may have been computed with higher precision and the result
may be there in higher precision, too. By telling the compiler
explicitly that I want the result converted to double (before
it is converted to int), I would expect the excess precision
to disappear.

So this is a workaround for a specific compiler/platform. It shouldn't make
a difference as far as the standard is concerned, right? Because if it did,
I'd be seriously confused.


Christian

 

[Michael Mair]

Michael Mair wrote:



So this is a workaround for a specific compiler/platform. It shouldn't make
a difference as far as the standard is concerned, right? Because if it did,
I'd be seriously confused.

Have a look at Lawrence Kirby's reply upthread:
"
How so? While gcc does have conformance issues in the area of floating
point I don't see any problem with the output of this program. In
particular C doesn't require n and nn to be set to the same value.
In fact the difference between them, where an intermediate value is
stored in a double object in one case and not in the other, is a classic
example of where the values are allowed to differ.
"
As Lawrence Kirby is nearly always right, I work from that assumption.
The cast IMO must give the same behaviour as using the intermediate
variable but as you can see in my reply, I phrased this as a question.


Cheers
Michael

 

[Kevin Bracey]

In message <[email protected]>

So this is a workaround for a specific compiler/platform. It shouldn't make
a difference as far as the standard is concerned, right? Because if it did,
I'd be seriously confused.

The expression "1.0 / d" has semantic type double but is allowed to contain
extra precision. The compiler is not required to round to double after every
calculation - this is to aid systems like the x86.

The macro FLT_EVAL_METHOD in <math.h> should be set to the correct value to
tell the programmer whether expressions are evaluated to extra precision in a
given implementation.

When you write n = 1.0 / d, the compiler is converting the high precision
result of the division directly to an int.

When you write n = (double) (1.0 / d), the compiler must convert the high
precision value to double precision, because the standard requires that extra
precision must be lost upon casting or assignment. That double is then
converted to an int, with potentially different results.

The compiler's actions in the OP's example are conforming. If it fails to
round correctly on a (double) cast, then that would be non-conforming. And
indeed, I believe many (all?) x86 compilers do fail to round excess precision
out in some circumstances. I wouldn't be surprised if the OP's program gave
different output if compiled with -O2 - the excess precision may fail to be
lost on assignment to dd.

 

[Christian Kandeler]

The expression "1.0 / d" has semantic type double but is allowed to
contain extra precision. The compiler is not required to round to double
after every calculation - this is to aid systems like the x86.

Wow, that's weird. I would certainly not expect an expression to change its
value when I cast it to the type it already has. Thanks for the
explanation.


Christian

 

[Mark McIntyre]

Wow, that's weird. I would certainly not expect an expression to change its
value when I cast it to the type it already has.

Its not changing its value - its discarding extra precision.

 

[Kevin Bracey]

In message <[email protected]>

Its not changing its value - its discarding extra precision.

Eh? How do you discard extra precision without changing the value? On such a
system,

1.0 / 3.0 == 1.0 / 3.0

would be true, but

1.0 / 3.0 == (double) (1.0 / 3.0)

would be false. That's changing the value, is it not? Or are you making some
subtle semantic point I'm missing?

 

[jacob navia]

good point!

 

[Randy Howard]

good point!

What is?

 

[jacob navia]

What is?

Kevin Bracey wrote:

Eh? How do you discard extra precision without changing the value? On such a
system,

1.0 / 3.0 == 1.0 / 3.0

would be true, but

1.0 / 3.0 == (double) (1.0 / 3.0)

would be false. That's changing the value, is it not? Or are you making some
subtle semantic point I'm missing?


That is a very good point. You can't discard the extra precision without
changing the value!

 

[Mark McIntyre]

That is a very good point. You can't discard the extra precision without
changing the value!

mathematically true, computationally not necessarily. Consider these
numbers

3.14159
3.1415926

Both, when stored in a 5-decimals register, would have the same value.

 

[Tim Rentsch]

Below is a program which converts a double to an integer in two
different ways, giving me two different values for the int. The basic
expression is 1.0 / (1.0 * 365.0) which should be 365, but one variable
becomes 364 and the other one becomes 365.

Does anyone have any insight to what the problem is?

Thanks in advance.
Bjørn

$ cat d.c
#include <stdio.h>

int main(void)
{
double dd, d = 1.0 / 365.0;
int n, nn;

n = 1.0 / d;
dd = 1.0 / d;
nn = dd;

printf("n==%d nn==%d dd==%f\n", n, nn, dd);
return 0;
}

$ gcc -Wall -O0 -ansi -pedantic -W -Werror -o d d.c
$ ./d
n==364 nn==365 dd==365.000000

$ gcc -v
Reading specs from /usr/lib/gcc/i386-redhat-linux/3.4.2/specs
Configured with: ../configure --prefix=/usr --mandir=/usr/share/man
--infodir=/usr/share/info --enable-shared --enable-threads=posix
--disable-checking --with-system-zlib --enable-__cxa_atexit
--disable-libunwind-exceptions --enable-java-awt=gtk
--host=i386-redhat-linux
Thread model: posix
gcc version 3.4.2 20041017 (Red Hat 3.4.2-6.fc3)
$

Having looked at what goes on with this sort of thing a lot during a
thread of several months ago, I feel obliged to post a response and
try to clear things up a bit.

A couple of postings have tried to sweep the problems under the rug
saying things like the problem is inherent in the nature of floating
point, or floating point isn't precise, or words to that effect. I
think that's both wrong and misleading. First, floating point is
*precise*; what it is not is *exact*. Floating point does behave
according to precise rules, which are implemention specific to some
extent, but they are precise nonetheless. The rules often don't give
the results we expect, because in mathematics we expect the results to
be exact, which floating point calculations are not. In spite of
that, floating point calculations do behave precisely, and if we
learn how they behave we can use them with confidence.

Furthermore, the question raised is not about the precision of the
answer but about whether floating point calculation is deterministic.
If you do the same calculation twice, do you get the same result? And
that question wasn't really addressed in other responses.

At least one posting attributed the problem to the x86 platform and
problems with gcc, with correcting/contradicting posts following. The
contradiction is both right and wrong. The behavior of the program
listed above is conformant with the C standard. But, this kind of
behavior does tend to show up erroneously on the x86 platform, and gcc
is known to have bugs around behavior of floating point (eg, 'double')
operands w.r.t. conforming to the C standard, especially in the
presence of optimization. It's probably worth bearing that in mind
when testing for how a program "should" behave. (See also the
workaround below.)

The key here, as was pointed out, is that the calculation '1.0 / d'
in the assignments

n = 1.0 / d;
dd = 1.0 / d;

is done (on the x86) in extended precision. In the second assignment,
the extended precision value is converted to 'double' before doing the
assignment; but in the first assignment, the extended precision value
is *not* converted to 'double' before being converted to 'int'. Thus,
in the (as was also posted) corrected code

n = (double) (1.0 / d);
dd = 1.0 / d;

we could reasonably expect that an assertion

assert( (int) dd == n );

to succeed. (Note: "reasonably expect" but not "certainly expect";
more below.)

It is the conversion to 'double' before the conversion to 'int' that
makes the two assigned expressions alike. The C standard requires
this; see 6.3.1.4, 6.3.1.5, and 6.3.1.8. In particular, either
casting with '(double)' or assigning to a 'double' variable is
required to convert an extended precision value to a value with
exactly 'double' precision. That's why the corrected code behaves
more as we expect.

That an assignment to a double variable sometimes requires a
conversion can lead to some unexpected results. For example, consider

#define A_EXPRESSION (*some side effect free expression*)
#define B_EXPRESSION (*some side effect free expression*)

double a = A_EXPRESSION;
double b = B_EXPRESSION;
double c = a + b;
double d = (A_EXPRESSION) + (B_EXPRESSION);

Normally we might expect that 'c' and 'd' can be used interchangeably
(when doing common subexpression elimination inside the compiler, for
example), but the conversion rules for C make this not so, or at least
not always so. Using temporaries in calculations using 'float' or
'double' changes the semantics in a way that doesn't happen with 'int'
values.

To return to the original question, how is it that the value of 'n'
differs from '(int) dd'? The behavior of the conversion to 'int' is
required by the standard to truncate. If the conversion from extended
precision to 'double' also truncated, there would be no discrepancy in
the program above. But the conversion (in this implementation) from
extended precision to 'double' doesn't truncate, it rounds. This
behavior conforms to the C standard, which requires that the result be
exact when possible, or either of the two nearest values otherwise
(assuming that there are such values represented in 'double'). Which
of the nearest values is chosen in the latter case is "implementation
defined".

So, as far as I can tell, an implementation could convert from
extended precision to 'double' by doing "statistical rounding" -- that
is, a mode in which rounding is done non-deterministically -- and
still be conformant. I don't know of any hardware that actually does
this, but as far as I can tell non-deterministic rounding is allowed.

Bottom line is, putting in the '(double)' cast will very likely make
the code deterministically yield 'n == (int) dd', but the standard
doesn't guarantee that it will. It may be possible to guarantee
deterministic behavior by setting the "rounding direction mode", or by
setting the "dynamic rounding precision mode", if supported; please
see Annex F.

======================================================================

Incidental note on getting around problems with gcc. Because gcc on
the x86 platform sometimes fails to convert extended precision values
to 'double' precision when the standard semantics require it, it's
nice to have a way to force a '(double)' conversion to take place.
The inline function

inline double
guarantee_double( double x ){
return * (volatile double *) &x;
}

is a pretty solid way of doing that.