float compares on x86 (i.e. 32 bit) bug?

[Matthias Kretz <kretz at fias dot uni-frankfurt dot de>]

Hi,

I was wondering why the following code failed:

int main() {
  const float_v x = ...;
  const float_v y = ...;
  float_v r = Vc::sqrt(x * x + y * y);

  for (int i = 0; i < Count; ++i) {
    const float a = r[0];
    const float b = std::sqrt(x[0] * x[0] + y[0] * y[0]);
    if (a != b) {
      std::cout << a << " != " << b << ": " << a - b << std::endl;
      return 1;
    }
  }
  return 0;
}

Output is:
1.41421 != 1.41421: 0

float_v is a class from my Vc library (http://gitorious.org/vc/) that wraps 
__m128. Vc::sqrt thus uses the SSE sqrt with single precision. std::sqrt OTOH 
is implemented as fsqrt and the comparison "a != b" is implemented with 
fucomip using the result from fsqrt directly and the result from Vc::sqrt 
loaded from memory. After fucomip the following flags are set:
(gdb) p $eflags
$6 = [ CF IF ]

Obviously I expected ZF to be set. Here's the problem:
st6 1.41421353816986083984375            (raw 0x3fffb504f30000000000)
st7 1.4142135623730950487637880730318329 (raw 0x3fffb504f333f9de6484)

Looking at the binary representation the single precision values are equal. 
Just st7 (the result from fsqrt) is not stored as single.

So now I understand why it failed. Am I right in my understanding that this is 
a bug? According to the C++ code the comparison should not have returned false 
because std::sqrt returns a single precision value.

Comments?

-- 
Matthias Kretz

Phone:   +49 6221 1871682
         sip:062211871682@sipgate.de
Email:   kretz@fias.uni-frankfurt.de