Trouble with integer floating point conversion

[rembremading]

Hi all!

The following piece of code has (for me) completely unexpected behaviour.
(I compile it with gcc-Version 4.0.3)
Something goes wrong with the integer to float conversion.
Maybe somebody out there understands what happens.
Essentially, when I subtract the (double) function value GRID_POINT(2) from
a variable which has been assigned the same value before this gives a
non-zero result and I really do not understand why.

The program prints
5.0000000000000 00000e-01; -2.7755575615628 91351e-17;
0.0000000000000 00000e+00

And the comparison
if(temp1==GRID_ POINT(2))
has negative outcome.

When I replace
((double)(i)) by 2.0
everything behaves as expected. So the output is
5.0000000000000 00000e-01; 0.0000000000000 00000e+00; 0.0000000000000 00000e+00

But: even if the integer to float conversion is inexact (which, I think,
should not be the case) something like
temp1 = GRID_POINT(2);
temp3 = temp1-GRID_POINT(2);
should still result in temp3==0.0, whatever function GRID_POINT does.

What do You think?

Thank you!

---------------------------------------------------
#include <stdio.h>

double GRID_POINT(int i);

double GRID_POINT(int i)
{
return ( 0.1 + ( (80.1-0.1)/(400.0) )*((double)(i)) );
}

int main (void) {

double temp1, temp2, temp3;

temp1 = GRID_POINT(2);
temp2 = GRID_POINT(2);
temp3 = temp1-GRID_POINT(2);

printf("%.18e; %.18e; %.18e\n", temp1, temp3, temp1-temp2 );

if(temp1==GRID_ POINT(2)){
printf("these two are equal\n");
}

return 0;
}
---------------------------------------------------

[Mark Bluemel]

Chris Torek wrote:

In article <fj**********@a ioe.org>
Mark Bluemel <ma**********@p obox.comwrote:

>... read http://www.network-theory.co.uk/docs...cintro_70.html

It answers your questions fairly comprehensively and gives a simple
way of disabling extended precision for the duration of a program's
execution.

It also notes, in passing, that setting the precision field in
the FPU control word does *not* affect the exponent range. This
means that the in-FPU computations are *still* not quite what
one might expect from an IEEE implementation.

I have the same note (with an example) at the
end of <http://web.torek.net/torek/c/numbers.html>.

As someone else noted else-thread, one can get "proper" IEEE floating
point on the x86 using gcc's "-ffloat-store" flag. This does,
however, slow down computation.

And still doesn't work for the OP's testcase :-

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

double GRID_POINT(int i) {
return ( 0.1 + ( (80.1-0.1)/(400.0) )*((double)(i)) );
}

int main (void) {
double temp1 = GRID_POINT(2);
double temp2 = temp1-GRID_POINT(2);

printf("%.18e; %.18e\n", temp1, temp2);

return 0;
}

$ cc -ffloat-store gr.c -o gr

$ ./gr
5.0000000000000 00000e-01; -2.7755575615628 91351e-17

As I commented else-thread -ffloat-store doesn't appear to alter how
a floating point return value from a function (presumably in a register)
is handled. This is probably covered by the comment in the manual page
:-

... "a few programs rely on the precise definition of IEEE floating
point. Use -ffloat-store for such programs, _after modifying them to
store all pertinent intermediate computations into variables_."

[Ben Bacarisse]

Mark Bluemel <ma**********@p obox.comwrites:

Chris Torek wrote:

<snip>

>As someone else noted else-thread, one can get "proper" IEEE floating
point on the x86 using gcc's "-ffloat-store" flag. This does,
however, slow down computation.

And still doesn't work for the OP's testcase :-

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

double GRID_POINT(int i) {
return ( 0.1 + ( (80.1-0.1)/(400.0) )*((double)(i)) );
}

int main (void) {
double temp1 = GRID_POINT(2);
double temp2 = temp1-GRID_POINT(2);

printf("%.18e; %.18e\n", temp1, temp2);

return 0;
}

$ cc -ffloat-store gr.c -o gr

$ ./gr
5.0000000000000 00000e-01; -2.7755575615628 91351e-17

A data point. It does here:

$ gcc -ffloat-store gr.c -o gr
$ ./gr
5.0000000000000 00000e-01; 0.0000000000000 00000e+00
$ gcc gr.c -o gr
$ ./gr
5.0000000000000 00000e-01; -2.7755575615628 91351e-17
$ gcc --version
gcc (GCC) 4.1.3 20070929 (prerelease) (Ubuntu 4.1.2-16ubuntu2)
<snip>

--
Ben.

[Alok Singhal]

On Thu, 13 Dec 2007 18:59:35 +0000, Mark Bluemel wrote:

Chris Torek wrote:

>>
As someone else noted else-thread, one can get "proper" IEEE floating
point on the x86 using gcc's "-ffloat-store" flag. This does, however,
slow down computation.

And still doesn't work for the OP's testcase :-

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

double GRID_POINT(int i) {
return ( 0.1 + ( (80.1-0.1)/(400.0) )*((double)(i)) );
}

int main (void) {
double temp1 = GRID_POINT(2);
double temp2 = temp1-GRID_POINT(2);

printf("%.18e; %.18e\n", temp1, temp2);

return 0;
}

$ cc -ffloat-store gr.c -o gr

$ ./gr
5.0000000000000 00000e-01; -2.7755575615628 91351e-17

This is OT, but the above happens because GRID_POINT(2) in the second call
is an 'unnamed temporary', and is not necessarily flushed to the memory by
gcc (even with -ffloat-store).

See http://arxiv.org/abs/cs/0701192 - "The pitfalls of verifying
floating-point computations".

-Alok

[Chris Torek]

>Chris Torek wrote:

>As someone else noted else-thread, one can get "proper" IEEE floating
point on the x86 using gcc's "-ffloat-store" flag. This does,
however, slow down computation.

In article <fj**********@a ioe.org>
Mark Bluemel <ma**********@p obox.comwrote:

>And still doesn't work for the OP's testcase ... [snip code; but it

reappears below, modified slightly]

>As I commented else-thread -ffloat-store doesn't appear to alter how
a floating point return value from a function (presumably in a register)
is handled.

This is allowed by the Standard -- you have to store the result
in a temporary variable, or use a cast, to trim off "extra"
precision.

>This is probably covered by the comment in the manual page :-
... "a few programs rely on the precise definition of IEEE floating
point. Use -ffloat-store for such programs, _after modifying them to
store all pertinent intermediate computations into variables_."

If gcc fails to store the number through memory (on the x86, that
is; this is how -ffloat-store trims the excess precision and
exponents) on a cast, that would be a bug.

On my Linux box here, the bug is not showing up, with the same
compilation options you used, or indeed any others. Then again,
I probably have a different Linux, different version of gcc, and
different libraries. (In particular I think this Linux and/or
libc defaults C programs to double, not extended, precision in
the FPU control word.)

Here is the modified test-case; compile with -DCAST_AWAY_EXCE SS to
test the effect of adding a cast.

% cat gr.c
#include <stdio.h>

double GRID_POINT(int i) {
return ( 0.1 + ( (80.1-0.1)/(400.0) )*((double)(i)) );
}

#ifdef USE_CAST
#define CAST_AWAY_EXCES S (double)
#else
#define CAST_AWAY_EXCES S /* nothing */
#endif

int main(void) {
double temp1 = GRID_POINT(2);
double temp2 = temp1 - CAST_AWAY_EXCES S GRID_POINT(2);

printf("%.18e; %.18e\n", temp1, temp2);

return 0;
}
% cc -o gr -ffloat-store gr.c
% ./gr
5.0000000000000 00000e-01; 0.0000000000000 00000e+00
--
In-Real-Life: Chris Torek, Wind River Systems
Salt Lake City, UT, USA (40°39.22'N, 111°50.29'W) [phone]
email: forget about it http://web.torek.net/torek/index.html
Reading email is like searching for food in the garbage, thanks to spammers.

[Mark Bluemel]

Ben Bacarisse wrote:

Mark Bluemel <ma**********@p obox.comwrites:

>Chris Torek wrote:

<snip>

>>As someone else noted else-thread, one can get "proper" IEEE floating
point on the x86 using gcc's "-ffloat-store" flag. This does,
however, slow down computation.

And still doesn't work for the OP's testcase :-

[Snip]

>
A data point. It does here:

$ gcc -ffloat-store gr.c -o gr
$ ./gr
5.0000000000000 00000e-01; 0.0000000000000 00000e+00
$ gcc gr.c -o gr
$ ./gr
5.0000000000000 00000e-01; -2.7755575615628 91351e-17
$ gcc --version
gcc (GCC) 4.1.3 20070929 (prerelease) (Ubuntu 4.1.2-16ubuntu2)

Looks like I need a later gcc then :-)

[CBFalconer]

James Kuyper wrote:

ym******@gmail. com wrote:

.... snip ...

>

>You quoted it: "How could the standard guarantee ANYTHING about
the precision of floating point calculations... "

Well, actually, the numerical limits section does specify the
precision. However, getting an inaccurate value with high
precision is not something most numerical analysts would be
happy about.

Please describe how you implement infinite precision in a finite
number of bits.

--
Merry Christmas, Happy Hanukah, Happy New Year
Joyeux Noel, Bonne Annee.
Chuck F (cbfalconer at maineline dot net)
<http://cbfalconer.home .att.net>

--
Posted via a free Usenet account from http://www.teranews.com

[Flash Gordon]

CBFalconer wrote, On 13/12/07 15:05:

James Kuyper wrote:

>ym******@gmail. com wrote:

... snip ...

>>You quoted it: "How could the standard guarantee ANYTHING about
the precision of floating point calculations... "

Well, actually, the numerical limits section does specify the
precision. However, getting an inaccurate value with high
precision is not something most numerical analysts would be
happy about.

Please describe how you implement infinite precision in a finite
number of bits.

That has nothing to do with James's comment. You do not need infinite
precision to get guaranteed accuracy. You do not need infinite precision
to get high precision.
--
Flash Gordon

[jameskuyper]

CBFalconer wrote:

James Kuyper wrote:

ym******@gmail. com wrote:

... snip ...

You quoted it: "How could the standard guarantee ANYTHING about
the precision of floating point calculations... "

Well, actually, the numerical limits section does specify the
precision. However, getting an inaccurate value with high
precision is not something most numerical analysts would be
happy about.

Please describe how you implement infinite precision in a finite
number of bits.

Do you have any particular reason for posing such an absurd challenge?
It's not relevant to any point that I, or anyone else, has made so far
in this discussion. Maybe if you explain more what you're looking for,
I might be able to help.

[Keith Thompson]

CBFalconer <cb********@yah oo.comwrites:

James Kuyper wrote:

>ym******@gmail. com wrote:

... snip ...

>>

>>You quoted it: "How could the standard guarantee ANYTHING about
the precision of floating point calculations... "

Well, actually, the numerical limits section does specify the
precision. However, getting an inaccurate value with high
precision is not something most numerical analysts would be
happy about.

Please describe how you implement infinite precision in a finite
number of bits.

I think you misunderstood James's use of the phrase "inaccurate
value".

For example, given a precision of 5 digits after the decimal point,
3.14159 is an accurate value of pi (to within the stated precision);
3.12345 is an innaccurate value of pi. Infinite precision is not
required for what's being discussed.

--
Keith Thompson (The_Other_Keit h) <ks***@mib.or g>
Looking for software development work in the San Diego area.
"We must do something. This is something. Therefore, we must do this."
-- Antony Jay and Jonathan Lynn, "Yes Minister"

[Barry Schwarz]

On Wed, 12 Dec 2007 16:09:12 +0100, rembremading
<re**********@g mx.netwrote:

>Hi all!

The following piece of code has (for me) completely unexpected behaviour.
(I compile it with gcc-Version 4.0.3)
Something goes wrong with the integer to float conversion.
Maybe somebody out there understands what happens.
Essentially, when I subtract the (double) function value GRID_POINT(2) from
a variable which has been assigned the same value before this gives a
non-zero result and I really do not understand why.

The program prints
5.000000000000 000000e-01; -2.7755575615628 91351e-17;
0.000000000000 000000e+00

And the comparison
if(temp1==GRID _POINT(2))
has negative outcome.

When I replace
((double)(i) ) by 2.0
everything behaves as expected. So the output is
5.000000000000 000000e-01; 0.0000000000000 00000e+00; 0.0000000000000 00000e+00

But: even if the integer to float conversion is inexact (which, I think,
should not be the case) something like
temp1 = GRID_POINT(2);
temp3 = temp1-GRID_POINT(2);
should still result in temp3==0.0, whatever function GRID_POINT does.

What do You think?

I think you should output the internal representation of the doubles
in hex so you can see the actual results of your computation, not what
printf thinks it should display.

>
Thank you!

---------------------------------------------------
#include <stdio.h>

double GRID_POINT(int i);

double GRID_POINT(int i)
{
return ( 0.1 + ( (80.1-0.1)/(400.0) )*((double)(i)) );

What purpose do you think the cast serves?

All the parentheses are superfluous except the pair around the
subtraction.

Do you really think 80.1 - .1 is exactly equal to 80.0?

Even if it is on your system, do you think 80.0/400.0 is exactly equal
to 0.2? Can any floating point value be exactly equal to 0.2 on your
system? (It is possible on mine but very few other people here have a
system with decimal floating point; most use binary.)

>}

int main (void) {

double temp1, temp2, temp3;

temp1 = GRID_POINT(2);
temp2 = GRID_POINT(2);
temp3 = temp1-GRID_POINT(2);

printf("%.18e; %.18e; %.18e\n", temp1, temp3, temp1-temp2 );

if(temp1==GRID_ POINT(2)){
printf("these two are equal\n");
}

return 0;
}
---------------------------------------------------

Remove del for email