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]

jacob navia wrote:

rembremading wrote:
[snip]

Can you tell me if you use the following program the result is better?

http://www.network-theory.co.uk/docs...cintro_70.html
(previously cited) gives a briefer version, which may be specific to
gcc, but that's what the OP has the issue with ....

[rembremading]

Yes, your assembler code works fine, although I don't understand what it
does :-)
For my purposes, however, it will be better to choose the constants such
that they are exactly representable, as outlined below.
Every additional function call, at this place in the code, is unwanted.
(although an inlined version of this function would be quite fast, I guess)
I would not want to activate the 64 bit computations for the entire program.
Actually for the most part of it calculation with higher precision are very
wellcome, as long as it does not affect the speed)
Nevertheless it is an interesting piece of code, which might be usefull in
the future.

jacob navia wrote:

rembremading wrote:
[snip]

Can you tell me if you use the following program the result is better?

This program sets the precision to 64 bits instead of
80 bits. (If I have no bugs... )

void SetPrecision64( void)
{
void (*fn)(void);
unsigned char code[] = {
0x50, // push %eax
0xd9,0x3c,0x24, // fnstcw (%esp,1)
0x8b,0x04,0x24, // mov (%esp,1),%eax
0x0f,0xba,0x34, 0x24,0x08, // btrl $0x8,(%esp,1)
0x66,0x81,0x0c, 0x24,0x00,0x02,// orw $0x200,(%esp,1)
0xd9,0x2c,0x24, // fldcw (%esp,1)
0x59, // pop %ecx
0xc3 // ret
};
fn = (void (*)(void))code;
fn();
}

[jacob navia]

rembremading wrote:

Yes, your assembler code works fine, although I don't understand what it
does :-)
For my purposes, however, it will be better to choose the constants such
that they are exactly representable, as outlined below.
Every additional function call, at this place in the code, is unwanted.
(although an inlined version of this function would be quite fast, I guess)
I would not want to activate the 64 bit computations for the entire program.
Actually for the most part of it calculation with higher precision are very
wellcome, as long as it does not affect the speed)
Nevertheless it is an interesting piece of code, which might be usefull in
the future.

You need it call it ONCE when your program STARTS.

Then, all calculations are done in 64 bits only
and not in 80.

That code sets the precision flag of the floating point unit
to 64 bits, i.e. double precision.

--
jacob navia
jacob at jacob point remcomp point fr
logiciels/informatique
http://www.cs.virginia.edu/~lcc-win32

[Mark Bluemel]

rembremading wrote:

Yes, your assembler code works fine, although I don't understand what it
does :-)

Then 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.

The mechanism outlined is the simplest way of "fixing" your problem.

[rembremading]

Yes, I understand that I would have to call it only once, but in the main
threefold for loop of my program I have to do something like

res = ( (1.0+f1)*(1.0+f 2)*f3*f4 - f1*f2*(1.0+f3)* (1.0+f4) );

Where f1..f4 are function values (not of the GRID_POINT function), which can
take values from 1.0e-100 to 1.0e2.
In this case, presumably, the intermediate values are stored in the register
and it is obvious that res is highly sensitive to the numerical precision
of f1..f4 and the intermediate results of the computation.
Therefore, if I can get the computation with 80 bit precision instead of 64
bit (at same costs), I definitely prefere the former one.
Although the final result is rounded to 64 bit precision this extra digits
can change a lot.

Because my system has dimension ~1e3 and the loop is threefold I cannot use
arbitrary precision math libraries. I cannot even use long double type
variables on my machine.

Now, the GRID_POINT function has to be evaluated when f1..f4 are computed.
Therefore I would have to switch 64 bit precision on (and off again) in
every single step. That is what I meant.

That code sets the precision flag of the floating point unit
to 64 bits, i.e. double precision.

I understood that, but I do not understand how the function does it
(possibly lack of assembler knowledge)

There might be also a lack of English knowledge, which could be the reason
for missunderstandi ngs. Sorry for that.

-Andreas
jacob navia wrote:

rembremading wrote:

>Yes, your assembler code works fine, although I don't understand what it
does :-)
For my purposes, however, it will be better to choose the constants such
that they are exactly representable, as outlined below.
Every additional function call, at this place in the code, is unwanted.
(although an inlined version of this function would be quite fast, I
guess) I would not want to activate the 64 bit computations for the
entire program. Actually for the most part of it calculation with higher
precision are very wellcome, as long as it does not affect the speed)
Nevertheless it is an interesting piece of code, which might be usefull
in the future.

You need it call it ONCE when your program STARTS.

Then, all calculations are done in 64 bits only
and not in 80.

That code sets the precision flag of the floating point unit
to 64 bits, i.e. double precision.

[jacob navia]

rembremading wrote:

Yes, I understand that I would have to call it only once, but in the main
threefold for loop of my program I have to do something like

res = ( (1.0+f1)*(1.0+f 2)*f3*f4 - f1*f2*(1.0+f3)* (1.0+f4) );

Where f1..f4 are function values (not of the GRID_POINT function), which can
take values from 1.0e-100 to 1.0e2.
In this case, presumably, the intermediate values are stored in the register
and it is obvious that res is highly sensitive to the numerical precision
of f1..f4 and the intermediate results of the computation.
Therefore, if I can get the computation with 80 bit precision instead of 64
bit (at same costs), I definitely prefere the former one.
Although the final result is rounded to 64 bit precision this extra digits
can change a lot.

Because my system has dimension ~1e3 and the loop is threefold I cannot use
arbitrary precision math libraries. I cannot even use long double type
variables on my machine.

Now, the GRID_POINT function has to be evaluated when f1..f4 are computed.
Therefore I would have to switch 64 bit precision on (and off again) in
every single step. That is what I meant.

>That code sets the precision flag of the floating point unit
to 64 bits, i.e. double precision.

I understood that, but I do not understand how the function does it
(possibly lack of assembler knowledge)

There might be also a lack of English knowledge, which could be the reason
for missunderstandi ngs. Sorry for that.

-Andreas

Your problem is a very complicated one Andreas. Maybe you will get
better advice
in sci.math.num-analysis, a group that is dedicated to this kind
of problems. I have seen discussions about problems of this type there,
and they have a math level that goes WAY beyond what we have here in
comp.lang.c.

I think you have to rearrange the terms in your equation to make them
less sensible to roundoff errors. I think in comp.math.num-analysis
you will find people that can help you to do that.
--
jacob navia
jacob at jacob point remcomp point fr
logiciels/informatique
http://www.cs.virginia.edu/~lcc-win32

[Keith Thompson]

ym******@gmail. com writes:

On Dec 13, 4:31 am, Keith Thompson <[e-mail] rote:

>jacob navia <ja...@nospam.c omwrites:

Flash Gordon wrote:

[...]

>>Since the standard doesn't force even IEEE754, there is *nothing*
the C language by itself can guarantee the user.

(I think jacob wrote the above; the attribution was snipped.)

>Are you seriously trying to claim that the only way to provide any
form of guarantee on floating point is to enforce IEEE754?

If there isn't even an accepted standard that is enforced, how
can you guarantee anything.

[...]

How could the standard guarantee ANYTHING about the precision of
floating point calculations when it doesn't even guarantee a common
standard?

Yes I am seriously saying that the absence of an enforced standard
makes any guarantee IMPOSSIBLE.

[...]

>Yes, jacob, it's true that the C standard makes no guarantees about
floating-point precision. It does make some guarantees about
floating-point range, which seems to contradict your claim above that
"there is *nothing* the C language by itself can guarantee the user".
(Perhaps in context it was sufficiently clear that you were talking
only about precision; I'm not going to bother to go back up the thread
to check.)

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

I was referring to jacob's earlier statement, which did not mention
precision:
| Since the standard doesn't force even IEEE754, there is *nothing*
| the C language by itself can guarantee the user.

>[...]

This is an example of how the "regulars" spread nonsense just with the
only objective of "contradict ing jacob" as they announced here
yesterday.

Tedious rant ignored.

So you just demonstrated that JN's paranoia is not *quite* paranoia.
Your post is exactly what he refers to (using strong words like
"nonsense", since JN gets pretty emotional in this grand war).

World will stop, and all newbies will get confused and die if you
stop arguing with Jacob, sure. You simlpy have to do it again and
again. Even when you don't have a point! Oh well.

Look again, I *agreed* with some of what jacob was saying; the
C standard doesn't make any guarantees about the precsion of
floating-point operations.

--
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"

[Keith Thompson]

Keith Thompson <ks***@mib.orgw rites:

ym******@gmail. com writes:

>On Dec 13, 4:31 am, Keith Thompson <[e-mail] rote:

[...]

>>Yes, jacob, it's true that the C standard makes no guarantees about
floating-point precision. It does make some guarantees about
floating-point range, which seems to contradict your claim above that
"there is *nothing* the C language by itself can guarantee the user".
(Perhaps in context it was sufficiently clear that you were talking
only about precision; I'm not going to bother to go back up the thread
to check.)

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

I was referring to jacob's earlier statement, which did not mention
precision:
| Since the standard doesn't force even IEEE754, there is *nothing*
| the C language by itself can guarantee the user.

And, to be clear, we should have been referring to accuracy rather
than precision.

--
Keith Thompson (The_Other_Keit h) ks***@mib.org <http://www.ghoti.net/~kst>
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"

[Keith Thompson]

Keith Thompson <ks***@mib.orgw rites:

ym******@gmail. com writes:

>On Dec 13, 4:31 am, Keith Thompson <[e-mail] rote:

[...]

>>Yes, jacob, it's true that the C standard makes no guarantees about
floating-point precision. It does make some guarantees about
floating-point range, which seems to contradict your claim above that
"there is *nothing* the C language by itself can guarantee the user".
(Perhaps in context it was sufficiently clear that you were talking
only about precision; I'm not going to bother to go back up the thread
to check.)

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

I was referring to jacob's earlier statement, which did not mention
precision:
| Since the standard doesn't force even IEEE754, there is *nothing*
| the C language by itself can guarantee the user.

And, to be clear, we should have been referring to accuracy rather
than precision.

--
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"

[Chris Torek]

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.
--
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.