18

I'm compiling and running the following program in 32 and 64 bit platforms:

int main()
{
  double y = 8.34214e08;
  double z = 1.25823e45;

  return y * z == 8.34214e08 * 1.25823e45;
}

While in 64bit the result is the expected (the values are the same and the exit code is non-zero) in 32bit seems there is a little difference between the value calculated at compile time, the right hand side of the comparison, and the left side computed at runtime.

Is this a bug in the compiler or there is a logical explanation?

EDIT: this is different from Why comparing double and float leads to unexpected result? because here all the values are double.

Community
  • 1
  • 1
waj
  • 1,185
  • 10
  • 14
  • 1
    @waj: I didn't downvote your question. It's a fine question. – Steve Summit Jul 14 '15 at 22:13
  • @SteveSummit ok, sorry about that then, I got the downvote at the same time you replied and StackOverflow doesn't show who voted. Anyway, I can live with it but I really need to know the reasons, that's why I'm sending a question here in the first place. – waj Jul 14 '15 at 22:16
  • @Giorgi: in that question there is a difference because the variable is declared as "float" and the literals are "double". Here all the values are double. – waj Jul 14 '15 at 22:19
  • @waj with gcc 5.1 it gives 1 on amd64 as well as armv7h (ARM 32 bit).. What compiler are you using? – woggioni Jul 14 '15 at 22:21
  • gcc appears to optimize the program down to `return 1` with any optimization level above `-O0`. – markgz Jul 14 '15 at 22:39
  • You have compared floats for equality. Go forth and sin no more. – Lee Daniel Crocker Jul 14 '15 at 22:59
  • @LeeDanielCrocker: :) I know, I'm not doing float equality on a real program. I just wanted to understand why the result is different when compiling and running on the same machine. – waj Jul 14 '15 at 23:06

3 Answers3

19

IEEE-754 allows intermediate computations to be done in a greater precision (emphasis mine).

(IEEE-754:2008) "A language standard should also define, and require implementations to provide, attributes that allow and disallow value-changing optimizations, separately or collectively, for a block. These optimizations might include, but are not limited to: [...] Use of wider intermediate results in expression evaluation."

In your case for example on a IA-32, the double values could be stored in the x87 FPU registers with greater precision (80-bit instead of 64). So you are actually comparing a multiplication done on double precision with a multiplication done on double-extended precision.

For example, on x64 where the result is 1 (the x87 FPU is not used as SSE is used instead), adding gcc option -mfpmath=387 to use the x87 makes the result change to 0 on my machine.

And if you wonder if that is also allowed by C, it is:

(C99, 6.3.1.p8) "The values of floating operands and of the results of floating expressions may be represented in greater precision and range than that required by the type;"

ouah
  • 142,963
  • 15
  • 272
  • 331
11

In general, never do equality checks with floating point numbers. You need to check whether the result you want differs from the result you get by less than a pre-set precision.

What is happening here is in all likelihood due to the multiplication being run on two different "platforms": once by your code, and once by the compiler, which may have a different precision. This happens with most compilers.

Your program would probably work if you compiled it with the same options that were used to compile the compiler (supposing the compiler was compiled by itself). But that would not mean you would get the correct result; you would be getting the same precision error the compiler is getting.

(Also, I'm assuming that the compiler performs a straight multiplication and the parsing code recognizing floats does not enter into the equation. This might well be wishful thinking on my part).

Testing

Using built-in specs.
COLLECT_GCC=gcc
COLLECT_LTO_WRAPPER=/usr/lib64/gcc/x86_64-suse-linux/4.8/lto-wrapper
Target: x86_64-suse-linux
Configured with: ../configure --prefix=/usr --infodir=/usr/share/info --mandir=/usr/share/man --libdir=/usr/lib64 --libexecdir=/usr/lib64 --enable-languages=c,c++,objc,fortran,obj-c++,java,ada --enable-checking=release --with-gxx-include-dir=/usr/include/c++/4.8 --enable-ssp --disable-libssp --disable-plugin --with-bugurl=http://bugs.opensuse.org/ --with-pkgversion='SUSE Linux' --disable-libgcj --disable-libmudflap --with-slibdir=/lib64 --with-system-zlib --enable-__cxa_atexit --enable-libstdcxx-allocator=new --disable-libstdcxx-pch --enable-version-specific-runtime-libs --enable-linker-build-id --enable-linux-futex --program-suffix=-4.8 --without-system-libunwind --with-arch-32=i586 --with-tune=generic --build=x86_64-suse-linux --host=x86_64-suse-linux
Thread model: posix
gcc version 4.8.3 20141208 [gcc-4_8-branch revision 218481] (SUSE Linux)



#include <stdio.h>

int main()
{
  double y = 8.34214e08;
  double z = 1.25823e45;

 return  printf("%s\n", y * z == 8.34214e08 * 1.25823e45 ? "Equal" : "NOT equal!");
}

Forcing -O0 to avoid the compiler from optimizing out the whole code (thanks @markgz!), we get

$ gcc -m32 -O0 -o float float.c && ./float NOT equal! $ gcc -m32 -frounding-math -O0 -o float float.c && ./float Equal

For the record, since you got there before me :-),

-frounding-math

Disable transformations and optimizations that assume default floating-point rounding behavior. This is round-to-zero for all floating point to integer conversions, and round-to-nearest for all other arithmetic truncations. This option should be specified for programs that change the FP rounding mode dynamically, or that may be executed with a non-default rounding mode. This option disables constant folding of floating-point expressions at compile time (which may be affected by rounding mode) and arithmetic transformations that are unsafe in the presence of sign-dependent rounding modes.

The default is -fno-rounding-math.

LSerni
  • 55,617
  • 10
  • 65
  • 107
  • Maybe I should clarify that I'm running the compilation and execution on the same machine each time. I totally agree that in general, equality should not be used for floating points. – waj Jul 14 '15 at 22:27
  • I was about to answer something like this, but what really puzzles me, is the difference between compile-time and run-time calculation. I would have thought GCC did the compile-time calculation using the same precision / settings as in the run-time calculation or whatever, but somehow the result differs :) I'm as puzzled as @waj – Morten Jensen Jul 14 '15 at 22:28
  • I remember getting weird results with `gcc` some years back, due to some math "optimization" I had overlooked. There are several "floating point" standards and gcc, for one, supports *all* of them. Of course gcc ran with only one - the one with which itself was compiled. And I might have been a bit, er, too free when building the toolchain. My professor was *really* not amused :-) – LSerni Jul 14 '15 at 22:38
  • You're right, I compiled again with -frounding-math and the results are now the same. I found more info here: https://gcc.gnu.org/wiki/FloatingPointMath#Note_on_built-in_math_functions – waj Jul 14 '15 at 22:42
  • 2
    Gcc fully optimizes the calculation down to 'true' with any optimization level above -O0 (due to constant folding). One would need to fool gcc's constant folding optimization with a more complex program to test the OP's original question. – markgz Jul 14 '15 at 23:00
5

Floating point calculations done at compile time often occur at a higher precision than double uses at run time. Also C may perform run-time intermediate double calculations at the higher long double precision. Either explain your inequality. See FLT_EVAL_METHOD for details.

  volatile double y = 8.34214e08;
  volatile double z = 1.25823e45;
  volatile double yz = 8.34214e08 * 1.25823e45;
  printf("%.20e\n", y);
  printf("%.20e\n", z);
  printf("%.20e\n", yz);
  printf("%.20Le\n", (long double) y*z);
  printf("%.20Le\n", (long double) 8.34214e08 * 1.25823e45);

8.34214000000000000000e+08
1.25822999999999992531e+45
// 3 different products!
1.04963308121999993395e+54
1.04963308121999993769e+54
1.04963308122000000000e+54

Your results may slightly differ.

chux - Reinstate Monica
  • 143,097
  • 13
  • 135
  • 256