GCC problem with raw double type comparisons

Question

I have the following bit of code, however when compiling it with GCC 4.4 with various optimization flags I get some unexpected results when its run.

#include <iostream>

int main()
{
   const unsigned int cnt = 10;
   double lst[cnt] = { 0.0 };
   const double v[4] = { 131.313, 737.373, 979.797, 731.137 };

   for(unsigned int i = 0; i < cnt; ++i) {
      lst[i] = v[i % 4] * i;
   }

   for(unsigned int i = 0; i < cnt; ++i) {
      double d = v[i % 4] * i;
      if(lst[i] != d) {
         std::cout << "error @ : " << i << std::endl;
         return 1;
      }
   }
   return 0;
}

• when compiled with: "g++ -pedantic -Wall -Werror -O1 -o test test.cpp" I get the following output: "error @ : 3"

• when compiled with: "g++ -pedantic -Wall -Werror -O2 -o test test.cpp" I get the following output: "error @ : 3"

• when compiled with: "g++ -pedantic -Wall -Werror -O3 -o test test.cpp" I get no errors

• when compiled with: "g++ -pedantic -Wall -Werror -o test test.cpp" I get no errors

I do not believe this to be an issue related to rounding, or epsilon difference in the comparison. I've tried this with Intel v10 and MSVC 9.0 and they all seem to work as expected. I believe this should be nothing more than a bitwise compare.

If I replace the if-statement with the following: if (static_cast<long long int>(lst[i]) != static_cast<long long int>(d)), and add "-Wno-long-long" I get no errors in any of the optimization modes when run.

If I add std::cout << d << std::endl; before the "return 1", I get no errors in any of the optimization modes when run.

Is this a bug in my code, or is there something wrong with GCC and the way it handles the double type?

Note: I've just tried this with gcc versions 4.3 and 3.3, the error is not exhibited.

Resolution: Mike Dinsdale noted the following bug report: http://gcc.gnu.org/bugzilla/show_bug.cgi?id=323 It seems the GCC team are not completely sure about nature of problem.

As suggested in the bug report a possible resolution is to use the ffloat-store option. I've tried this and it works, however the results from a performance point of view are not that great, though ymmv.

Answer 1

The problem is likely the result of losing some precision when storing the result of an expression vs. the compiler not doing so in a local as an optimization:

double d = v[i % 4] * i;  // the result, `d`, might be kept in a register 
                          //   instead of storing it in a memory location, 
                          //   keeping full precision

if(lst[i] != d) {         // the value stored in lst[i] may have lost some 
                          //   precision since it had to be stored in memory,
                          //   which might not be able to hold the full 
                          //   precision that the expression generated

The C99 standard says in 6.3.1.8/2 "Usual arithmetic conversions":

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; the types are not changed thereby.

Answer 2

The fact that the result depends on the optimization settings suggests it might be the x87 extended precision messing with things (as Michael Burr says).

Here's some code I use (with gcc on x86 processors) to switch the extended precision off:

static const unsigned int PRECISION_BIT_MASK = 0x300;
///< bitmask to mask out all non-precision bits in the fpu control word \cite{INTEL}.
static const unsigned int EXTENDED_PRECISION_BITS = 0x300;
///< add to the fpu control word (after zeroing precision bits) to turn on extended precision \cite{INTEL}.
static const unsigned int STANDARD_PRECISION_BITS = 0x200;
///< add to the fpu control word (after zeroing precision bits) to turn off extended precision \cite{INTEL}.

void set_fpu_control_word(unsigned int mode)
{
  asm ("fldcw %0" : : "m" (*&mode));
}

unsigned int get_fpu_control_word()
{
  volatile unsigned int mode = 0;
  asm ("fstcw %0" : "=m" (*&mode));
  return mode;
}

bool fpu_set_extended_precision_is_on(bool state)
{
  unsigned int old_cw = get_fpu_control_word();
  unsigned int masked = old_cw & ~PRECISION_BIT_MASK;
  unsigned int new_cw;
  if(state)
    new_cw = masked + EXTENDED_PRECISION_BITS;
  else
    new_cw = masked + STANDARD_PRECISION_BITS;
  set_fpu_control_word(new_cw);
  return true;
}

bool fpu_get_extended_precision_is_on()
{
  unsigned int old_cw = get_fpu_control_word();
  return  ((old_cw & PRECISION_BIT_MASK) == 0x300);
}

Or you can just run your code with valgrind, which doesn't simulate the 80-bit registers, and is probably easier for a short program like this!

Answer 3

The width of the floating point registers in x86 is different from the width of the double in RAM. Therefore comparisons may succeed or fail depending entirely on how the compiler decides to optimize the loads of floating point values.