Is this GCC optimization incorrect?

Question

I have a functor which takes a value, casts it to double, takes the log and casts the value back to the original type. For the purpose of this question, the original and output type is float. Here is the original C++ code:

return static_cast< TOutput >( std::log( static_cast< double >( A ) ) ) 

When I compile in debug mode, everything goes as expected and GCC calls the underlying log function:

  51:/myfile.h ****     return static_cast< TOutput >( std::log( static_cast< double >( A ) ) );  219133                 .loc 112 51 0  219134 0010 488B45F0     movq  -16(%rbp), %rax # A, tmp64  219135 0014 F30F1000     movss (%rax), %xmm0 # *A_1(D), D.237346  219136 0018 0F14C0       unpcklps  %xmm0, %xmm0  # D.237346, D.237346  219137 001b 0F5AC0       cvtps2pd  %xmm0, %xmm0  # D.237346, D.237347  219138 001e E8000000     call  log #  219138      00  219139 0023 660F14C0     unpcklpd  %xmm0, %xmm0  # D.237347  219140 0027 660F5AC0     cvtpd2ps  %xmm0, %xmm0  # D.237347, D.237346  219141 002b F30F1145     movss %xmm0, -20(%rbp)  # D.237346, %sfp  219141      EC  219142 0030 8B45EC       movl  -20(%rbp), %eax # %sfp, <retval> 

However, when I turn optimizations on (-O2 -ggdb3 -DNDEBUG), it calls the logf (???) function:

  51:/myfile.h ****     return static_cast< TOutput >( std::log( static_cast< double >( A ) ) );  145171                 .loc 64 51 0  145172 01a0 F30F1004     movss (%rdx,%rax,4), %xmm0  # MEM[(const float &)_84], MEM[(const float &)_84]  145172      82           145173 01a5 E8000000     call  logf  # 

It gives a different output. Is this normal? Am I doing anything wrong? It seems to me that GCC is taking a very liberal interpretation of my code, which I wouldn't expect in the absence of the -ffast-math option.

Answer 1

It is a borderline optimization to transform the conversion to float of an application of the double-precision log to a float to an application of the single-precision log, but it can be argued to be acceptable.

Assuming that logf is correctly rounded and that the double-precision log is also correctly rounded or at least faithfully-rounded, the two computations will rarely differ. They can differ (for some rare inputs) because of double-rounding (in which “double” means “twice” and does not refer to the type). Double-rounding is statistically all the less significant that there are additional digits in the intermediate type's significand compared to the final type's significand (and this statistical argument is slightly rubbish from a mathematical point of view, but it “works in practice” for functions that have not been designed to be counter-examples). For didactic reasons, people (or Wikipedia) explain it with one or two extra digits of precision, but when you have 53 - 24 = 29 extra binary digits, it can be expected to happen as rarely as once in 2²⁹.

I am surprised by the optimization, and I would disturbed by it if I had written the code myself for an exhaustive search of double-rounding problems with log, but considering that the C++ standard does not mandate any level of accuracy from std::log, it is possible to consider it “not a bug”.

---

If, instead of log, we had been talking of one of the basic operations (such as *), then the transformation would have been incorrect, for a compiler that claims to provide IEEE 754 semantics, when it introduces visible changes. For a basic operation, the accuracy is specified indirectly by IEEE 754, and the specification leaves no room for variations.

It so happens that for the basic operations, there cannot be a visible change when floatOP(flx,fly) replaces (float)doubleOP((double)flx, (double)fly) (this thesis demonstrates this in chapter 6) , but there can be visible differences when the types are double and long double. This exact bug was recently fixed in Clang by Stephen Canon.

Answer 2

Yes, this optimisation is incorrect. log and logf are only required to be faithfully-rounded, so one could have

logf(4) = 0x1.62e42p+0 log(4)  = 0x1.62e42fefa39efp+0 

Changing the upconversion, log, and downconversion to a logf call may give incorrect results.

Pascal Cuoq's answer correctly points out that, if logf is correctly-rounded and log isn't garbage, then the results probably won't differ. However, logf on my platform is not correctly-rounded:

logf(0x1.306p-138)       = -0x1.7decc8p+6 (float)log(0x1.306p-138) = -0x1.7decc6p+6 mpfr_log(0x1.306p-138)   = -0x1.7decc6ff8a7a4a4450e9p+6 

Thankfully, I'm unable to reproduce this "optimisation" with my gcc.