5

The discussion started under my answer to another question. The following code determines machine epsilon:

float compute_eps() {
  float eps = 1.0f;

  while (1.0f + eps != 1.0f)
    eps /= 2.0f;

  return eps;
}

In the comments it was proposed that the 1.0f + eps != 1.0f test might fail because C++ standard permits the use of extra precision. Although I'm aware that floating-point operations are actually performed in higher precision (than specified by the actual types used), I happen to disagree with this proposal.

I doubt that during the comparison operations, such as == or !=, the operands are not truncated to the precision of their type. In other words, 1.0f + eps can of course be evaluated with the precision higher than float (for example, long double), and the result will be stored in the register that can accommodate long double. However, I think that before performing the != operation left operand will be truncated from long double to float, hence the code can never fail to determine eps precisely (i.e. it can never do more iterations than intended).

I haven't found any clue on this particular case in C++ standard. Furthermore, the code works fine and I'm sure that the extra precision technique is used during its execution because I have no doubt that any modern desktop implementation in fact uses extra precision during calculations.

What do you think about it?

Community
  • 1
  • 1
Alexander Shukaev
  • 16,674
  • 8
  • 70
  • 85
  • Not truncating to the precision of the type before comparison is exactly what some GCC users are complaining about in several duplicates of bug 323: http://gcc.gnu.org/bugzilla/show_bug.cgi?id=323 (I don't remember the specific bug reports, sorry, bug 323 is a bit of a grab-all). One example is like `double d = 3. / 7.; if (d != 3. / 7.) /* taken */ …` – Pascal Cuoq May 01 '13 at 19:52
  • Why can't you just use `std::numeric_limits::epsilon()`? – masrtis May 01 '13 at 19:54
  • @masrtis: That's off-topic. I can, the point is different here. – Alexander Shukaev May 01 '13 at 19:56
  • I haven't looked for it (it may be from the C standard), but my recollection is that the rule is that floating-point calculations can be done at higher precision, and that **storing** the resulting value into a variable of the appropriate type reduces the precision to that of the stored type. And even then, some compilers persist with higher precision unless you use a command-line option that says to respect the rules on floating-point stores. Either way, a comparison is not a store, so either or both sides of the comparison can be calculated at higher precision. – Pete Becker May 01 '13 at 20:00

1 Answers1

4

Sorry that this example is C and not C++. It should not be difficult to adapt:

~ $ gcc -mfpmath=387 -mno-sse2  c.c
~ $ ./a.out 
incredible but true.
~ $ gcc -mfpmath=sse -msse2  c.c
~ $ ./a.out 
~ $ cat c.c
#include "stdio.h"

double d = 3. / 7.;
double d1 = 3.;

int main() {
  if (d != d1 / 7.)
    printf("incredible but true.\n");
  return 0;
}

gcc -msse2 -mfpmath=sse is a strict IEEE 754 compiler. With that compiler, the if is never taken. However, gcc -mno-sse2 -mfpmath=387 has to use the 387 unit with its higher precision. It does not reduce the precision before the != test. The test ends up comparing the extended-precision result of 3. / 7. on the right-hand side to the double-precision result of the same division on the left-hand side. This cause a behavior that may appear strange.

Both gcc -msse2 -mfpmath=sse and gcc -mno-sse2 -mfpmath=387 are standard-compliant. It is only the case that the former has it easy, generating SSE2 instructions, and thus can provide a strict IEEE 754 implementation, whereas the latter has to do its best with an ancient instruction set.

A loop such as:

while (eps1 != 1.0f)
  eps /= 2.0f, eps1 = 1.0f + eps;

with eps1 declared of type float should be more robust with respect to extended precision.

The compiler that generates x87 code that does not truncate before comparison is this one:

~ $ gcc -v
Using built-in specs.
Target: i686-apple-darwin11
Configured with: /private/var/tmp/llvmgcc42/llvmgcc42-2336.11~148/src/configure --disable-checking --enable-werror --prefix=/Applications/Xcode.app/Contents/Developer/usr/llvm-gcc-4.2 --mandir=/share/man --enable-languages=c,objc,c++,obj-c++ --program-prefix=llvm- --program-transform-name=/^[cg][^.-]*$/s/$/-4.2/ --with-slibdir=/usr/lib --build=i686-apple-darwin11 --enable-llvm=/private/var/tmp/llvmgcc42/llvmgcc42-2336.11~148/dst-llvmCore/Developer/usr/local --program-prefix=i686-apple-darwin11- --host=x86_64-apple-darwin11 --target=i686-apple-darwin11 --with-gxx-include-dir=/usr/include/c++/4.2.1
Thread model: posix
gcc version 4.2.1 (Based on Apple Inc. build 5658) (LLVM build 2336.11.00)

Here is another:

~ $ clang -mno-sse2  c.c
~ $ ./a.out 
incredible but true.
~ $ clang -v
Apple LLVM version 4.2 (clang-425.0.24) (based on LLVM 3.2svn)
Target: x86_64-apple-darwin12.3.0
Thread model: posix
Pascal Cuoq
  • 79,187
  • 7
  • 161
  • 281
  • Well, does that mean that my assumptions are correct? I.e. if everything is standard compliant, and no subtle extensions or bugs are present, then the operands are truncated before the comparison? – Alexander Shukaev May 01 '13 at 20:07
  • @Haroogan If the operands were truncated before the comparison, the message “incredible but true.” would never show. It shows for some compilers, including one whose characteristics I have given in my post. – Pascal Cuoq May 01 '13 at 20:09
  • What if we add explicit cast: `float(1.0f + eps) != 1.0f`? Does that amend the issue? – Alexander Shukaev May 01 '13 at 21:31
  • @Haroogan Yes, a cast should also remove excess precision. See Joseph S. Myers's analysis of how an implementation should conform: http://gcc.gnu.org/ml/gcc-patches/2008-11/msg00105.html – Pascal Cuoq May 01 '13 at 21:55