Precidents, Precision, or Plain Screwed Up?

[buddyletts]

Below is some test code:

#include <iostream>
#include <cmath>
using namespace std;

int main()
{
double P = -1943.6650940829325;
double R = 3.6683591763449312;
double Rc = 3.5964170000000002;
double a = exp(-pow(R,2)/pow(Rc,2));
double b = a*P;
double c = P*a;
double d = P*exp(-pow(R,2)/pow(Rc,2));
double e = exp(-pow(R,2)/pow(Rc,2))*P;
P *= exp(-pow(R,2)/pow(Rc,2));

return 0;
}

I'm looking for the value of P. (Above code written for VC++ 6.0.
Library modifications required for the HP.)

On a HP c3750 running aCC: P = -686.71739393609221
On an iMac G5 running Tiger w/gcc: P = -686.71739393609221
On a P3 running ME/ MS VC++ 6.0: P = -686.71739393609209

On the same P3 w/MS VC++ 6.0: b = -686.71739393609221
c = -686.71739393609221
d = -686.71739393609209
e = -686.71739393609209

So my question is why are b and c different than d, e, and P if the
exponent is evaluated before the multiplication? Any comments would be
appreciated. Thanks in advance.

 

[Victor Bazarov]

Below is some test code:

#include <iostream>
#include <cmath>
using namespace std;

int main()
{
double P = -1943.6650940829325;
double R = 3.6683591763449312;
double Rc = 3.5964170000000002;
double a = exp(-pow(R,2)/pow(Rc,2));
double b = a*P;
double c = P*a;
double d = P*exp(-pow(R,2)/pow(Rc,2));
double e = exp(-pow(R,2)/pow(Rc,2))*P;
P *= exp(-pow(R,2)/pow(Rc,2));

return 0;
}

I'm looking for the value of P. (Above code written for VC++ 6.0.
Library modifications required for the HP.)

On a HP c3750 running aCC: P = -686.71739393609221
On an iMac G5 running Tiger w/gcc: P = -686.71739393609221
On a P3 running ME/ MS VC++ 6.0: P = -686.71739393609209

On the same P3 w/MS VC++ 6.0: b = -686.71739393609221
c = -686.71739393609221
d = -686.71739393609209
e = -686.71739393609209

BTW, what's the "correct" answer?

Also BTW, I ran your example on P4 after compiling it with VC6 and VC7.1
and got ...09 _in_all_vars_ from VC6 and ...21 _in_all_vars_ from VC7.1.

So my question is why are b and c different than d, e, and P if the
exponent is evaluated before the multiplication?

They weren't when I ran it.

> Any comments would be
appreciated. Thanks in advance.

<implementation-specific>

If the exponent is evaluated and then _stored_ in a 'double' variable,
you *may* lose several significant digits. All calculations in the x87
are done in a 10-byte "extended double precision" format. If there is
no intervening store operation, the extra digits at the end of the value
are kept around. That's what happens when 'd' and 'e' are calculated,
I believe. Try using the compile switch "improve FP consistency".

</implementation-specific>

I suggest you post your question to microsoft.public.vc.language for
more information on compiler switches and FP behaviour.

V

 

[red floyd]

Victor Bazarov wrote:
[redacted]

<implementation-specific>

If the exponent is evaluated and then _stored_ in a 'double' variable,
you *may* lose several significant digits. All calculations in the x87
are done in a 10-byte "extended double precision" format. If there is
no intervening store operation, the extra digits at the end of the value
are kept around. That's what happens when 'd' and 'e' are calculated,
I believe. Try using the compile switch "improve FP consistency".

</implementation-specific>

I suggest you post your question to microsoft.public.vc.language for
more information on compiler switches and FP behaviour.

V

<implementation-specific>
Also, AFAICT, VC (through 7.1) makes its "long double" 64 bits, the same
as a regular double. There's no way that I can tell to get VC7.1 to use
an 80 bit native double.
</implementation-specific>

 

[Lionel B]

Below is some test code:

#include <iostream>
#include <cmath>
using namespace std;

int main()
{
double P = -1943.6650940829325;
double R = 3.6683591763449312;
double Rc = 3.5964170000000002;
double a = exp(-pow(R,2)/pow(Rc,2));
double b = a*P;
double c = P*a;
double d = P*exp(-pow(R,2)/pow(Rc,2));
double e = exp(-pow(R,2)/pow(Rc,2))*P;
P *= exp(-pow(R,2)/pow(Rc,2));

return 0;
}

I'm looking for the value of P. (Above code written for VC++ 6.0.
Library modifications required for the HP.)

On a HP c3750 running aCC: P = -686.71739393609221
On an iMac G5 running Tiger w/gcc: P = -686.71739393609221
On a P3 running ME/ MS VC++ 6.0: P = -686.71739393609209

On the same P3 w/MS VC++ 6.0: b = -686.71739393609221
c = -686.71739393609221
d = -686.71739393609209
e = -686.71739393609209

So my question is why are b and c different than d, e, and P if the
exponent is evaluated before the multiplication? Any comments would
be appreciated. Thanks in advance.

Floating point is not an exact representation - it can only store real numbers with limited precision. Thus floating
point calculation on a computer is *not* equivalent to mathematical calculation with real numbers. This may have
dismaying (to the uninitiated) consequences for results of floating point arithmetic. In particular:

(i) quantities calculated via different methods by floating point arithmetic which you think should be identical on
mathematical grounds will not necessarily be so in practice

(ii) the same calculation on different systems will not necessarily deliver identical results.

and, BTW:

(iii) you have no guarantee that library functions like exp(), pow(), etc. even use the same algorithm, let alone
deliver results to the same precision on different systems.

Have a good read of the David Goldberg article "What Every Computer Scientist Should Know About Floating-Point
Arithmetic":

http://docs.sun.com/source/806-3568/ncg_goldberg.html

 

[buddyletts]

I realize the above and have evaluated exp and pow functions on all
machines as well as all of the terms and expressions in between. As
far as I can tell I get identical results on all three platforms for
every instance of each of the terms and expressions EXCEPT for when I
put everything on 1 line on the P3 running VC++ 6.0 (i.e. the
expressions: d, e, and P).

BTW, when evaluated on the Mac and HP all term b, c, d, e, & P are
identical. It's only the terms d, e, & P running VC++ 6.0 on the P3
that I get something different.

Also BTW, the Mac, HP, and b&c terms on the P3 are the more correct
answer (-686.71739393609221). When evaluating the expression using
variable precision correct answer is -686.717393936092225...

I guess my question is more along the lines of: What's the difference
between storing the exponent in a variable and then multiplying it vs.
just flat out doing the multiplication on one line? (And that's
probably not the best way to word it.) On the HP and Mac the answer
appears to be nothing. On the P3 running VC++ 6.0, the answer appears
to be 0.00000000000012.

 

[Lionel B]

[snip]

I guess my question is more along the lines of: What's the difference
between storing the exponent in a variable and then multiplying it vs.
just flat out doing the multiplication on one line?

I suspect the answer may lie in the storing of intermediates with unpredictable precision - see Victor Basarov's reply
earlier in this thread. Have you tried playing around with your compilers' switches for enforcing IEEE fp compliance?
Eg. the -ffloat-store switch on GCC:

Quote:

"Do not store floating point variables in registers, and inhibit other options that might change whether a floating
point value is taken from a register or memory.
This option prevents undesirable excess precision on machines such as the 68000 where the floating registers (of the
68881) keep more precision than a double is supposed to have. Similarly for the x86 architecture. For most programs, the
excess precision does only good, but 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."

HTH,