CODEWITHSUNDEEP
c++linuxgccassemblyfloating-point
static_rtti
static_rtti

Reputation: 56262

Is this GCC optimization incorrect?

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.

Upvotes: 24

Views: 2184

Answers (4)

Pascal Cuoq
Pascal Cuoq

Reputation: 80276

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 229.

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.

Upvotes: 33

tmyklebu
tmyklebu

Reputation: 14205

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.

Upvotes: 2

SylvainL
SylvainL

Reputation: 3938

You shouldn't have any concern about the different output. If you cast a float value to a double and then back to a float value; you can get a different value; therefore the different value that you get from casting a float to a double value, performing the log operation and then casting down the result back to a float is not better than the value that you get by calling directly the logf function.

If you have any concern about the precision of the operation, then you shouldn't use the float type at all and all your mathematical values should be double all the way.

By using the float type, you signal the compilator that you don't have any special attention to the final precision of your mathematical operations beyond that of the one provided by the float type. Put in other words, when you do some mathematical operations with real numbers (numbers other than integers), the precision can only either remains the same or decrease; it can never increase; therefore, if you start with a float type, you cannot increase your precision beyond that type.

Upvotes: -1

user3185968
user3185968

Reputation:

I've tried an equivalent C program:

#include <math.h>
#include <stdlib.h>
#include <stdio.h>


static double dostuff(float in)
{
  return log(in);
}

int main(int argc, char **argv)
{
  if (argc < 2)
    exit(EXIT_FAILURE);
  float in = strtof(argv[1], NULL);
  float out = dostuff(in);
  printf("%f\n", out);
  return 0;
}

And found that gcc does not use logf even when using -Ofast. The only thing -Ofast enables is __log_finite() to be used. However, changing the return type of dostuff() to float does enable this optimization.

Upvotes: 1

Related Questions