By using this site, you agree to our updated Privacy Policy and our Terms of Use. Manage your Cookies Settings.
424,680 Members | 1,868 Online
Bytes IT Community
+ Ask a Question
Need help? Post your question and get tips & solutions from a community of 424,680 IT Pros & Developers. It's quick & easy.

double to int conversion yields strange results

P: n/a
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.
Bjrn

$ 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)
$
Nov 14 '05 #1
Share this Question
Share on Google+
31 Replies


P: n/a
=?ISO-8859-1?Q?Bj=F8rn_Augestad?= <bo*@metasystems.no> wrote:
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.


How about people reading the FAQ for a change?

<http://www.eskimo.com/~scs/C-faq/q14.1.html>
<http://www.eskimo.com/~scs/C-faq/q14.4.html>
<http://www.eskimo.com/~scs/C-faq/q14.5.html>

Richard
Nov 14 '05 #2

P: n/a
<http://www.eskimo.com/~scs/C-faq/q14.5.html>


Even though it says, a==b is wrong, some _simple_ like

double i = strtod(argv[1], NULL);
if(i == 1.0) {
//enjoy
}

seems to work. Is it due to the nature that 1.0 can be encoded "good enough"
in FP regs?

Jan Engelhardt
--
If you knew the language, you'd know that: 私は考える、従って私はある。
Nov 14 '05 #3

P: n/a
Richard Bos wrote:
=?ISO-8859-1?Q?Bj=F8rn_Augestad?= <bo*@metasystems.no> wrote:

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.

How about people reading the FAQ for a change?

<http://www.eskimo.com/~scs/C-faq/q14.1.html>
<http://www.eskimo.com/~scs/C-faq/q14.4.html>
<http://www.eskimo.com/~scs/C-faq/q14.5.html>

Richard


Thanks, but why do you think this applies to my question?

Bjrn
Nov 14 '05 #4

P: n/a


Bjrn Augestad wrote:
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.
Bjrn

$ 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)
$


gcc on x86 architectures may keep certain double variables in
registers where they are effectively long doubles and some may
be stored in memory where they really are doubles (talking about
64 bits vs 80 bits).
Try the option
-ffloat-store
and, if that does not "help", also
-mfpmath=sse -msse2
or whatever is appropriate for you.
This may be inefficient, so you could also use the macro
#define EVAL_AS_DOUBLE(d) (*((volatile double *)&(d)))
on your double objects, whenever you need to be sure that double
really is double.
Cheers
Michael
--
E-Mail: Mine is a gmx dot de address.

Nov 14 '05 #5

P: n/a
Michael Mair wrote:


Bjrn Augestad wrote:
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?
[snip]

gcc on x86 architectures may keep certain double variables in
registers where they are effectively long doubles and some may
be stored in memory where they really are doubles (talking about
64 bits vs 80 bits).
Try the option
-ffloat-store
and, if that does not "help", also
-mfpmath=sse -msse2
or whatever is appropriate for you.
This may be inefficient, so you could also use the macro
#define EVAL_AS_DOUBLE(d) (*((volatile double *)&(d)))
on your double objects, whenever you need to be sure that double
really is double.


Thanks, Michael.
-msse2 did the trick,-ffloat-store and -mfpmath didn't work.

What an odd bug that was. I had been using -march=prescott and
everything was working fine. Then I removed that option for portability
reasons and then the code started to break in strange ways.

I'll rewrite the code to avoid the problems completely.
Bjrn
Nov 14 '05 #6

P: n/a
=?ISO-8859-1?Q?Bj=F8rn_Augestad?= <bo*@metasystems.no> wrote:
Richard Bos wrote:
=?ISO-8859-1?Q?Bj=F8rn_Augestad?= <bo*@metasystems.no> wrote:
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.


How about people reading the FAQ for a change?

<http://www.eskimo.com/~scs/C-faq/q14.1.html>
<http://www.eskimo.com/~scs/C-faq/q14.4.html>
<http://www.eskimo.com/~scs/C-faq/q14.5.html>


Thanks, but why do you think this applies to my question?


_Never_ assume that floating point computations are exact. You are
obviously computing 364.999999999 in one case, and 365.000000000001 in
the other. If you'd read the FAQ, you'd have realised this.

Richard
Nov 14 '05 #7

P: n/a

Richard Bos wrote:
=?ISO-8859-1?Q?Bj=F8rn_Augestad?= <bo*@metasystems.no> wrote:

Richard Bos wrote:
=?ISO-8859-1?Q?Bj=F8rn_Augestad?= <bo*@metasystems.no> wrote:
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.

How about people reading the FAQ for a change?

<http://www.eskimo.com/~scs/C-faq/q14.1.html>
<http://www.eskimo.com/~scs/C-faq/q14.4.html>
<http://www.eskimo.com/~scs/C-faq/q14.5.html>


Thanks, but why do you think this applies to my question?

_Never_ assume that floating point computations are exact. You are
obviously computing 364.999999999 in one case, and 365.000000000001 in
the other. If you'd read the FAQ, you'd have realised this.


However, this is a gcc bug due to the incredibly *** x86 architecture.
Bjorn got the semantics right, so even though pointing people to the
FAQ is hardly ever wrong, his problem was of a different nature and
could not be solved within the program because gcc left the realms of
standard C.
Cheers
Michael
--
E-Mail: Mine is a gmx dot de address.

Nov 14 '05 #8

P: n/a
Bjrn Augestad wrote:
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?
That you, like many clueless people, have not learned how floating point
works, have not bothered to check the FAQ before posting, have not
checked the archives for countless threads on the same subject -- some
in the last week -- started by other clueless wonders, or even followed
the newsgroup for a week before posting. Will September *never* end?


--
"If you want to post a followup via groups.google.com, don't use
the broken "Reply" link at the bottom of the article. Click on
"show options" at the top of the article, then click on the
"Reply" at the bottom of the article headers." - Keith Thompson
Nov 14 '05 #9

P: n/a
Bjrn Augestad wrote:

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?

$ 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;
}


After all the huffing and puffing about the FAQ, nobody has told
him what the fundamental cause is. The x86 does floating point
calculations with an 80 bit operand, which allows for extra
precision. Those operands are in the FP processor internal stack,
and stay there until needed. When the computation 1.0 / d is
directly converted to an integer, the extra bits are present. When
it is first stored in a 64 bit double, they have to be removed, and
thus the final conversion to an integer is different.

There are various options to control this, but fundamentally the
FAQ is correct - floating point values are imprecise. Optimization
level will probably also affect the result.

--
"If you want to post a followup via groups.google.com, don't use
the broken "Reply" link at the bottom of the article. Click on
"show options" at the top of the article, then click on the
"Reply" at the bottom of the article headers." - Keith Thompson
Nov 14 '05 #10

P: n/a
On Fri, 11 Feb 2005 13:58:12 +0100, Michael Mair wrote:

....
However, this is a gcc bug due to the incredibly *** x86 architecture.
Bjorn got the semantics right, so even though pointing people to the
FAQ is hardly ever wrong, his problem was of a different nature and
could not be solved within the program because gcc left the realms of
standard C.


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.

Lawrence
Nov 14 '05 #11

P: n/a
This happens because gcc sets the default precision to full
precision in the x86, i.e. 80 bits precision.

Using the lcc-win32 compiler with your program I obtain the
same result as gcc.

If I add however a call to a function to set the precision
of all operations to 64 bits only, I get the correct
result:
#include <stdio.h>
#include <fenv.h>
int main(void)
{
double dd, d = 1.0 / 365.0;
int n, nn;

DoublePrecision(); // Set the machine to 64 bits
// precision *only*
n = (double)(1.0 / d);
dd = 1.0 / d;
nn = dd;

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

n==365 nn==365 dd==365.000000

If you do not use the lcc-win32 compiler, gcc will
swallow this assembly:
.text
.globl _DoublePrecision
_DoublePrecision:
push %eax
fnstcw (%esp)
movl (%esp),%eax
btr $8,(%esp)
orw $0x200,(%esp)
fldcw (%esp)
popl %ecx
movb %ah,%al
andl $3,%eax
ret

You compile it with gcc and you set the machine precision
at the beginning of the prograM.

Maybe gcc offers a functional equivalent, I do not know.
Nov 14 '05 #12

P: n/a


Lawrence Kirby wrote:
On Fri, 11 Feb 2005 13:58:12 +0100, Michael Mair wrote:

...

However, this is a gcc bug due to the incredibly *** x86 architecture.
Bjorn got the semantics right, so even though pointing people to the
FAQ is hardly ever wrong, his problem was of a different nature and
could not be solved within the program because gcc left the realms of
standard C.


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.


Sorry, I shot off the other mail thinking of the following situation:

Q: If the code had been
double dd, d = 1.0 / 365.0;
int n, nn;

n = (double) (1.0 / d); /*difference*/
dd = 1.0 / d;
nn = dd;
I would expect n==nn. Is there a loop hole in the standard
s.th. this could be not the case?

TIA
Michael
--
E-Mail: Mine is a gmx dot de address.

Nov 14 '05 #13

P: n/a
On Fri, 11 Feb 2005 12:46:30 +0100, in comp.lang.c , Jan Engelhardt
<je*****@linux01.gwdg.de> wrote:
<http://www.eskimo.com/~scs/C-faq/q14.5.html>


Even though it says, a==b is wrong, some _simple_ like

if(i == 1.0) {

seems to work. Is it due to the nature that 1.0 can be encoded "good enough"
in FP regs?


Possibly, tho there's no guarantee. It might work today, and fail
tomorrow., or work on one platform but not on another. Or always work
perfectly on every platform you're able to test on *except* the one your
new customer is going to use but hasn't told you about yet....

--
Mark McIntyre
CLC FAQ <http://www.eskimo.com/~scs/C-faq/top.html>
CLC readme: <http://www.ungerhu.com/jxh/clc.welcome.txt>
Nov 14 '05 #14

P: n/a
Michael Mair wrote:
...
_Never_ assume that floating point computations are exact. You are
obviously computing 364.999999999 in one case, and 365.000000000001 in
the other. If you'd read the FAQ, you'd have realised this.


However, this is a gcc bug due to the incredibly *** x86 architecture.
...


This has absolutely nothing to do with x86 architecture. These issues
are caused by certain properties of the floating-point representation
(natural to virtually any architecture) and by some inherent properties
of binary notation system itself. For example, floating-point binary
notation of a decimal value '0.1' has fractional part of infinite
length, which means that this value cannot be represented precisely by
any architecture using traditional floating-point representation.

In OP's case, 1 is divided by 365. The expected resultant value is
1/365, which is

0.(000000001011001110001100111110011011)

in binary notation. It is a periodic fraction of infinite length. It
cannot be ever represented precisely in floating-point format of finite
length on any architecture. Ever. Absolutely impossible. And it has
nothing to do with x86 and/or GCC.

There are alternative ways to represent fractional numbers, which can at
least partially solve these issues, by they have their own drawbacks.

--
Best regards,
Andrey Tarasevich
Nov 14 '05 #15

P: n/a
In article <42***********************@news.wanadoo.fr>
jacob navia <ja***@jacob.remcomp.fr> wrote:
This happens because gcc sets the default precision to full
precision in the x86, i.e. 80 bits precision.
Actually, it does not: it leaves the precision bits alone. Someone
or something else is responsible for setting them initially.
If I add however a call to a function to set the precision
of all operations to 64 bits only, I get the correct result ...
"Correct" according to some other standard than Standard C,
anyway. (Standard C is pretty lax about FP results.)

Unfortunately, setting the precision control word on the x86 is
not sufficient in general. It will solve this particular problem,
but not other problems on the x86.
If you do not use the lcc-win32 compiler, gcc will
swallow this assembly:
(I assume you mean gas -- gcc itself is a C compiler, not an
assembler. :-) )

[snippage]
Maybe gcc offers a functional equivalent, I do not know.


GCC has a way to do the trick with inline assembly, but this is
off-topic, and does not solve the problem anyway, because the x86
FPU is obnoxious. (Exercise for the x86 student: what 64-bit result
do you expect for 1.0e200 * 1.0e200? What 80-bit result do you
expect for 1.0e200L * 1.0e200L? What 80-bit result will you get
in the FPU stack for the first expression, no matter how you fiddle
with the control word? Does it matter if overflows do not occur
when they should?)
--
In-Real-Life: Chris Torek, Wind River Systems
Salt Lake City, UT, USA (4039.22'N, 11150.29'W) +1 801 277 2603
email: forget about it http://web.torek.net/torek/index.html
Reading email is like searching for food in the garbage, thanks to spammers.
Nov 14 '05 #16

P: n/a
Andrey Tarasevich wrote:
Michael Mair wrote:
...
_Never_ assume that floating point computations are exact. You are
obviously computing 364.999999999 in one case, and 365.000000000001 in
the other. If you'd read the FAQ, you'd have realised this.


However, this is a gcc bug due to the incredibly *** x86 architecture.
...

This has absolutely nothing to do with x86 architecture. These issues
are caused by certain properties of the floating-point representation
(natural to virtually any architecture) and by some inherent properties
of binary notation system itself. For example, floating-point binary
notation of a decimal value '0.1' has fractional part of infinite
length, which means that this value cannot be represented precisely by
any architecture using traditional floating-point representation.

In OP's case, 1 is divided by 365. The expected resultant value is
1/365, which is

0.(000000001011001110001100111110011011)

in binary notation. It is a periodic fraction of infinite length. It
cannot be ever represented precisely in floating-point format of finite
length on any architecture. Ever. Absolutely impossible. And it has
nothing to do with x86 and/or GCC.

There are alternative ways to represent fractional numbers, which can at
least partially solve these issues, by they have their own drawbacks.


Please have a look at my reply to Lawrence Kirby's reply.
I was thinking about the wrong situation and thought the OP
had used a cast -- in which case I would have expected the same
result but would not have gotten it with gcc.

-Michael
--
E-Mail: Mine is an /at/ gmx /dot/ de address.
Nov 14 '05 #17

P: n/a
Chris Torek wrote:

GCC has a way to do the trick with inline assembly, but this is
off-topic, and does not solve the problem anyway, because the x86
FPU is obnoxious. (Exercise for the x86 student: what 64-bit result
do you expect for 1.0e200 * 1.0e200? What 80-bit result do you
expect for 1.0e200L * 1.0e200L? What 80-bit result will you get
in the FPU stack for the first expression, no matter how you fiddle
with the control word? Does it matter if overflows do not occur
when they should?)


Overflow will be detected only when WRITING the numbers...

This precision setting will be used when writing/reading numbers
into the FPU, not when doing operations, as it should...

That is why the precision setting doesn't provoke an
overflow when multiplying 1e200*1e200: the number fits
in 80 bits with a dynamic range until around 1.18e4932,
only when you attempt to write the result into RAM
does the CPU notice... Ohh shh... I was in overflow!

jacob
Nov 14 '05 #18

P: n/a
Michael Mair wrote:
n = (double) (1.0 / d); /*difference*/


How does this differ from "n*=*1.0*/*d"?
Christian
Nov 14 '05 #19

P: n/a


Christian Kandeler wrote:
Michael Mair wrote:

n = (double) (1.0 / d); /*difference*/


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.
In fact, if I try to calculate ***_EPSILON with gcc
and test
fltepsilon = 1.0;
while(1.0F < 1.0+fltepsilon)
fltepsilon /= 2.0;
fltepsilon *= 2.0;
I find fltepsilon==LDBL_EPSILON but if I insert a cast
while(1.0F < (float)(1.0+fltepsilon))
I arrive at the right result. The same does not hold for double
and DBL_EPSILON unless I force gcc to do so.

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.
Cheers
Michael
--
E-Mail: Mine is a gmx dot de address.

Nov 14 '05 #20

P: n/a
Michael Mair wrote:


Christian Kandeler wrote:
Michael Mair wrote:

n = (double) (1.0 / d); /*difference*/

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.

Bjrn

Nov 14 '05 #21

P: n/a
Michael Mair wrote:
n = (double) (1.0 / d); /*difference*/


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.


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
Nov 14 '05 #22

P: n/a


Christian Kandeler wrote:
Michael Mair wrote:

n = (double) (1.0 / d); /*difference*/

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.


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
--
E-Mail: Mine is a gmx dot de address.

Nov 14 '05 #23

P: n/a
In message <37*************@individual.net>
Christian Kandeler <ch****************@hob.de_invalid> wrote:
Michael Mair wrote:
n = (double) (1.0 / d); /*difference*/

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.


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.

--
Kevin Bracey, Principal Software Engineer
Tematic Ltd Tel: +44 (0) 1223 503464
182-190 Newmarket Road Fax: +44 (0) 1728 727430
Cambridge, CB5 8HE, United Kingdom WWW: http://www.tematic.com/
Nov 14 '05 #24

P: n/a
Kevin Bracey wrote:
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
Nov 14 '05 #25

P: n/a
On Mon, 14 Feb 2005 16:13:20 +0100, in comp.lang.c , Christian Kandeler
<ch****************@hob.de_invalid> wrote:
Kevin Bracey wrote:
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.


Its not changing its value - its discarding extra precision.

--
Mark McIntyre
CLC FAQ <http://www.eskimo.com/~scs/C-faq/top.html>
CLC readme: <http://www.ungerhu.com/jxh/clc.welcome.txt>

----== Posted via Newsfeeds.Com - Unlimited-Uncensored-Secure Usenet News==----
http://www.newsfeeds.com The #1 Newsgroup Service in the World! 120,000+ Newsgroups
----= East and West-Coast Server Farms - Total Privacy via Encryption =----
Nov 14 '05 #26

P: n/a
In message <km********************************@4ax.com>
Mark McIntyre <ma**********@spamcop.net> wrote:
On Mon, 14 Feb 2005 16:13:20 +0100, in comp.lang.c , Christian Kandeler
<ch****************@hob.de_invalid> wrote:
Kevin Bracey wrote:
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.


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?

--
Kevin Bracey, Principal Software Engineer
Tematic Ltd Tel: +44 (0) 1223 503464
182-190 Newmarket Road Fax: +44 (0) 1728 727430
Cambridge, CB5 8HE, United Kingdom WWW: http://www.tematic.com/
Nov 14 '05 #27

P: n/a
good point!
Nov 14 '05 #28

P: n/a
In article <42***********************@news.wanadoo.fr>,
ja***@jacob.remcomp.fr says...
good point!


What is?

--
Randy Howard (2reply remove FOOBAR)
"Making it hard to do stupid things often makes it hard
to do smart ones too." -- Andrew Koenig
Nov 14 '05 #29

P: n/a
Randy Howard wrote:
In article <42***********************@news.wanadoo.fr>,
ja***@jacob.remcomp.fr says...
good point!

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!
Nov 14 '05 #30

P: n/a
On Wed, 16 Feb 2005 11:15:09 +0100, in comp.lang.c , jacob navia
<ja***@jacob.remcomp.fr> wrote:
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.
--
Mark McIntyre
CLC FAQ <http://www.eskimo.com/~scs/C-faq/top.html>
CLC readme: <http://www.ungerhu.com/jxh/clc.welcome.txt>
Nov 14 '05 #31

P: n/a
Bjrn Augestad <bo*@metasystems.no> writes:
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.
Bjrn

$ 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.
Nov 14 '05 #32

This discussion thread is closed

Replies have been disabled for this discussion.