C++ double precision failure with Visual Studio 2012 / Intel Compiler under Windows

Question

i have a problem with the double precision under Windows with the Visual Studio and the Intel Compiler.

#include <stdio.h>
#include <math.h>
int main()
{
  double result = 42;
  int rounds = 14;
  for(int a=1; a<=rounds; a++)
  {
       result = sqrt(sqrt(result));
       printf("Round %2.i is %.17g\n", a, result);
   }
   for(int a=1; a<=rounds; a++)
   {
       result = pow(result, 4);
       printf("Round %2.i is %.17g\n", a, result);
   }
  return 0;
}

With this code i produce under Linux with these compilersettings:

Intel Compiler under linux:

icc -O2 -fp-model double test.cpp -o test

GCC Compiler under linux:

g++ -O2 test.cpp -o test

and Matlab and Scilab under Windows the correct result:

Round  1 is 2.5457298950218306
Round  2 is 1.2631446315995756
Round  3 is 1.0601401197016873
Round  4 is 1.0147073765340913
Round  5 is 1.0036567375971333
Round  6 is 1.0009129334669549
Round  7 is 1.0002281552726189
Round  8 is 1.0000570339386641
Round  9 is 1.0000142581797196
Round 10 is 1.0000035645258711
Round 11 is 1.0000008911302767
Round 12 is 1.0000002227824947
Round 13 is 1.000000055695619
Round 14 is 1.0000000139239045
Round  1 is 1.000000055695619
Round  2 is 1.0000002227824947
Round  3 is 1.0000008911302765
Round  4 is 1.0000035645258705
Round  5 is 1.0000142581797171
Round  6 is 1.0000570339386543
Round  7 is 1.0002281552725802
Round  8 is 1.0009129334668005
Round  9 is 1.0036567375965142
Round 10 is 1.0147073765315875
Round 11 is 1.0601401196912235
Round 12 is 1.263144631549705
Round 13 is 2.545729894619797
Round 14 is 41.999999973468661

With this Intel compileroptions under linux

icc -O2 test.cpp -o test

and this gcc compileroptions under linux

g++ -O2 -ffast-math test.cpp -o test

and all Visual Studio 2012 Compiler Settings i get this wrong result:

Round  1 is 2.5457298950218306
Round  2 is 1.2631446315995756
Round  3 is 1.0601401197016873
Round  4 is 1.0147073765340913
Round  5 is 1.0036567375971333
Round  6 is 1.0009129334669549
Round  7 is 1.0002281552726189
Round  8 is 1.0000570339386641
Round  9 is 1.0000142581797196
Round 10 is 1.0000035645258711
Round 11 is 1.0000008911302767
Round 12 is 1.0000002227824947
Round 13 is 1.000000055695619
Round 14 is 1.0000000139239045
Round  1 is 1.000000055695619
Round  2 is 1.0000002227824947
Round  3 is 1.0000008911302767
Round  4 is 1.0000035645258711
Round  5 is 1.0000142581797198
Round  6 is 1.0000570339386647
Round  7 is 1.0002281552726222
Round  8 is 1.0009129334669684
Round  9 is 1.0036567375971872
Round 10 is 1.0147073765343093
Round 11 is 1.0601401197025984
Round 12 is 1.2631446316039172
Round 13 is 2.5457298950568319
Round 14 is 42.000000002309839

Under Windows the following settings are ignored and i always get the wrong result.

Intel compileroptions under windows:

icl /O2 /fp:double test.cpp -o test.exe

Visual Studio compileroptions under Windows:

/fp:precise
/fp:strict
/fp:fast

produced all the same "wrong numbers"...

Why is this error and how can i fix it that i get on all platforms the same working double precision?

Thank you very much and greetings

Natalie

Edit:

To benchmark i have looped the code:

#include <stdio.h>
#include <math.h>
int main()
{
  double x = 1;
  double result = 0;
  int rounds = 12;
  while (x!=result)
  {
  x++;
  result=x;
  for(int a=1; a<=rounds; a++)
  {
       result = sqrt(sqrt(result));
   }
   for(int a=1; a<=rounds; a++)
   {
       result = pow(result, 4);
   }
   }
   printf("The next possible with %2.i rounds is %.0f\n", rounds, result);
  return 0;
}

Compiled with the working linux compilers: Intel Compiler:

icc -O2 -fp-model double speedtest3.cpp -o speedtest3

Result:

The next possible with 12 rounds is 3671078
real    0m1.950s
user    0m1.947s
sys     0m0.003s

GCC Compiler:

g++ -O2 speedtest3.cpp -o speedtest3

Result:

The next possible with 12 rounds is 3671078

real    0m3.445s
user    0m3.442s
sys     0m0.004s

Code with pow:

#include <stdio.h>
#include <math.h>
int main()
{
  double x = 1;
  double result = 0;
  int rounds = 12;
  while (x!=result)
  {
  x++;
  result=x;
  for(int a=1; a<=rounds; a++)
  {
       result = pow(result, 0.25);
   }
   for(int a=1; a<=rounds; a++)
   {
       result = pow(result, 4);
   }
   }
   printf("The next possible with %2.i rounds is %.0f\n", rounds, result);
  return 0;
}

Compiled with the same options: Result Intel Compiler:

The next possible with 12 rounds is 3671078

real    0m2.887s
user    0m2.885s
sys     0m0.004s

Result GCC Compiler:

The next possible with 12 rounds is 3671078

real    0m5.905s
user    0m5.905s
sys     0m0.008s

Results:

sqrt is near 2x faster then pow The intel compiler is near 2x faster then the gcc compiler

Answer 1

They are both correct.

When dealing with floating point math at any precision, you cannot demand that the result be any exact number to the last digit. Expecting to get exactly 41.999999973468661 is just as misguided as expecting to get exactly 42.000000000000000.

Here's an example of how you might get different results: double-precision floating point arithmetic in C++ is (typically, depends on the compiler and platform) 64 bits. However, Intel CPU's floating point hardware calculates numbers internally using 80 bits. At some point, the values are pulled from the FPU and trimmed into 64-bits. What "at some point" means depends on what the compiler and it's optimizer thinks is appropriate or optimal, and therefore, different compilers/settings will result in rounding at different parts of the computation, affecting the exact digits of the result.

Even using the same exact binary program on different computers might legally change the result, and still be "correct". For example, it used to be that some intel processors didn't have a hardware floating point unit at all (ages ago, the last one being the 486SX). In those cases, your program would run the computations on a separate software-only library that would do nothing at all in 80-bit precision.

Addendum:

I want to put things in perspective. Wolfram Alpha says the distance between New York and San Francisco is about 2577 Miles. At that scale, the discrepancy you are so worried about between the two numbers is like saying that your own estimate of the distance "is wrong by 1/10^th inch" ^(*). If that divergence is significant to your problem, floating point arithmetic is the wrong tool for the job.

---

^* Relative error is (42.000000002309839 - 41.999999973468661)/41.999999973468661 = 6.8669471471949833966026905617771e-10. 2577 Miles = 163278720 inches. The equivalent divergence on that measurement would be (163278720 inches * 6.8669471471949833966026905617771e-10) = 0.1121226... inches. All calculations done in Windows Calculator on Windows 7, which uses fixed-point arithmetic.

Answer 2

Look at one way pow(double,double) is implemented, here. Different compilers from different vendors are bound to contain different implementations, with resulting differences.

As has been suggested, a better option is to use

double h = x*x;
x = h*h; 

What you call "incorrect" is more accurate, compare:

delta is 2.3098394308362913e-09   42.00000000230983 - 42
delta is 2.6531338903623691e-08   42 - 41.999999973468661

Later What are you trying to prove or show? Executing some numeric computation repeatedly is usually worse than computing the result directly; therefore pow( 42, 1.0/(1024*1024*64)) is theoretically preferable to doing taking 4th roots in a loop, and so is the reverse process: pow( result, 1024*1024*64 ), which might enable pow to use the minimum number of multiplications.

Anyway, the umpteenth root of 42 cannot be represented as a machine number, and so the umpteenth power of that is bound to deviate from 42. You can estimate the error using elementary numerical math.