C/C++ and Fortran

Originally posted by: Jinsong JI

---

Hey, when I am porting my code from GCC to XLC, I got some floating point precision error.

I spent days to narrow down the issue and seems like there is a potential xlc bug when doing floating point calculation! 
Can you please have a look?
Attached the simple testcase and reproducing steps:

$cat maf.c

#include <stdio.h>
double a, b, c;
int main()
{
  a = 1.60618e-07;
  b = -0.000128960;
  c = 0.00124548;
  c = a + (b * c);
  printf("%.16E\n",c);
  return (0);
}

$cat t1.ksh

echo "xlc"
xlc maf.c -o xl "$@" && ./xl
echo "gcc"
gcc maf.c -o g  && ./g

$./t1.ksh

xlc
8.9919999997675296E-13
gcc
8.9919999996900178E-13

 

---

 

 It is not the first time we receive bug report like this. The testcase did show the fact that XLC is producing different result from GCC.

 Is this really a xlc bug?

No! This is not a bug, it is actually caused by different floating point support strategy between different compilers.

By default, XL C/C++ follows most but not all of the rules in the IEEE standard.

 As we konw, floating point calculation result may vary with different rounding mode.

 In this case, whether we do rounding after b*c will definitely affect the final result, and it is the root cause of the above testcase!

 As we can see in http://www-01.ibm.com/support/knowledgecenter/SSXVZZ_13.1.2/com.ibm.xlcpp1312.lelinux.doc/proguide/maf.html

Handling multiply-add operations

By default, the compiler generates a single non-IEEE 754 compatible multiply-add instruction for binary floating-point expressions such as a+b*c, partly because one instruction is faster than two. Because no rounding occurs between the multiply and add operations, thismight also produce a more precise result. However, the increased precision might lead to different results from those obtained in other environments, and might cause x*y-x*y to produce a nonzero result. To avoid these issues, you can suppress the generation of multiply-add instructions by using the -qfloat=nomaf option."

 Let us see the test result with -qfloat=nomaf.
 

$ ./t1.ksh -qfloat=nomaf

xlc
8.9919999996900178E-13
gcc
8.9919999996900178E-13

  Now, we can get the same precision with GCC ! (although it is less precise now).
  

  There might be some other small difference with floating point in default setting, especially when you use -O3 and above.  

  xlc/xlf may do some transformation that are essential to delivering greater performance at higher optimization levels. [1]

  For example:

• optimization of x + 0.0 to x (because x might be -0.0, and this would change the sign)
• reordering floating point expressions
• rewriting x / 4.5 as x * (1/4.5) (because the lower order bit of the result may differ)
• replacing a loop of calls to sin with a call to vector_sin
• use of VMX vector float instructions 

  Most programs will give the same results, but it does matter for a few programs, and for that reason,

  there is a power option -qstrict that can help you to do trade-off between precision and performance.

  We would recommend you read through "Handling floating-point operations" in Optimization and Programming Guide before you start porting floating point applications.

   http://www-01.ibm.com/support/knowledgecenter/SSXVZZ_13.1.2/com.ibm.xlcpp1312.lelinux.doc/proguide/cvriscfp.html