Bug 57974 - std::pow(std::complex<long double>(0),1) returns (-nan,-nan)
Summary: std::pow(std::complex<long double>(0),1) returns (-nan,-nan)
Status: UNCONFIRMED
Alias: None
Product: gcc
Classification: Unclassified
Component: middle-end (show other bugs)
Version: 4.9.0
: P2 normal
Target Milestone: ---
Assignee: Not yet assigned to anyone
URL:
Keywords: wrong-code
Depends on:
Blocks:
 
Reported: 2013-07-24 21:56 UTC by Henner Sudek
Modified: 2013-08-18 12:11 UTC (History)
3 users (show)

See Also:
Host:
Target:
Build:
Known to work:
Known to fail:
Last reconfirmed:


Attachments
the preprocessed file (*.i*) that triggers the bug (87.51 KB, text/plain)
2013-07-24 21:56 UTC, Henner Sudek
Details

Note You need to log in before you can comment on or make changes to this bug.
Description Henner Sudek 2013-07-24 21:56:58 UTC
Created attachment 30548 [details]
the preprocessed file (*.i*) that triggers the bug

The following program gives the wrong result when compiled with the options given below. Removing any of the compiler options gives the expected result (0,0). I was able to reproduce the error with gcc-4.7.3, 4.8.1 and 4.9.0.

$ cat complex.cpp
#include <complex>
#include <iostream>

int main(){
  std::cerr << std::pow(std::complex<long double>(0),1) << std::endl;
}

$ c++ -std=c++11 -fno-math-errno -funsafe-math-optimizations complex.cpp
$ ./a.out 
(-nan,-nan)

$ /home/henner/gcc-4.9.0/bin/g++ -v 
Using built-in specs.
COLLECT_GCC=/home/henner/gcc-4.9.0/bin/g++
COLLECT_LTO_WRAPPER=/home/henner/gcc-4.9.0/libexec/gcc/x86_64-unknown-linux-gnu/4.9.0/lto-wrapper
Target: x86_64-unknown-linux-gnu
Configured with: ../configure --prefix=/home/henner/gcc-4.9.0
Thread model: posix
gcc version 4.9.0 20130724 (experimental) (GCC)
Comment 1 Paolo Carlini 2013-07-24 23:05:08 UTC
IMHO the most interesting bit is that it only happens for long double. Anyway, this is neither C++ library nor front-end, either middle-end or libc, likely the former because on my machine neither clang++ nor icc has the issue while the underlying libc is the same glibc. Anyway, reduced:

  __builtin_expl (- 1 / 0.0)
Comment 2 Paolo Carlini 2013-07-24 23:17:04 UTC
Or, more elegantly:

  __builtin_expl (-__builtin_huge_vall ())
Comment 3 Paolo Carlini 2013-07-25 15:22:11 UTC
Uros, can you help me with this annoying wrong-code?

I suspect the back end could be also involved because it happens only for long double and I didn't see anything special for long double in builtins.c

(if isn't clear enough, the expression in Comment #1 or #2 evaluates to -nan instead of 0)
Comment 4 Uroš Bizjak 2013-07-25 15:42:11 UTC
(In reply to Paolo Carlini from comment #3)

> I suspect the back end could be also involved because it happens only for
> long double and I didn't see anything special for long double in builtins.c

-funsafe-math-optimizations are used here:

'-funsafe-math-optimizations'

     Allow optimizations for floating-point arithmetic that (a) assume
     that arguments and results are valid and (b) may violate IEEE or
     ANSI standards.  When used at link-time, it may include libraries
     or startup files that change the default FPU control word or other
     similar optimizations.

The __builtin_expl is expanded to:

        flds    .LC0(%rip)
        fld     %st(0)
        frndint
(*)     fsubr   %st, %st(1)
        fxch    %st(1)
        f2xm1
        fadds   .LC2(%rip)
        fscale
        fstp    %st(1)
        fstpt   16(%rsp)

and with fsubr, we hit Inf - Inf, which is by definition NaN.

Obviously, this testcase exposes the unsafe part of -funsafe-math-optimizations. Without this option everything works OK.
Comment 5 Paolo Carlini 2013-07-25 15:47:40 UTC
But isn't this a bug? I mean, naively, what do we gain from the optimization point of view from not evaluating as 0 in any case? And why it happens only for long double?
Comment 6 Uroš Bizjak 2013-07-25 15:58:13 UTC
(In reply to Paolo Carlini from comment #5)
> But isn't this a bug? I mean, naively, what do we gain from the optimization
> point of view from not evaluating as 0 in any case? And why it happens only
> for long double?

Because long double always goes through x87. doubles and floats will fail your test with -m32 or -mfpmath=387 in addition to -funsafe-math-optimizations.

BTW: The test does use infinty (which is IEEE thingy), and there is explicit warning that -funsafe-math-optimizations will violate these rules.
Comment 7 Paolo Carlini 2013-07-25 16:08:55 UTC
However it's still not clear to me why this inconsistency doesn't happen with clang or icc, for example. I'm not convinced we are doing our job in the best way and I don't think we are going to make users happy. It's pretty easy to internally call in the c++ library or in user code another builtin which does produce a -inf without the user knowing and then passes it to a function like exp. If some of our builtins produce something which our other builtins aren't able to handle consistently we are looking for trouble.
Comment 8 Uroš Bizjak 2013-07-25 16:12:35 UTC
(In reply to Paolo Carlini from comment #7)
> However it's still not clear to me why this inconsistency doesn't happen
> with clang or icc, for example. I'm not convinced we are doing our job in
> the best way and I don't think we are going to make users happy. It's pretty
> easy to internally call in the c++ library or in user code another builtin
> which does produce a -inf without the user knowing and then passes it to a
> function like exp. If some of our builtins produce something which our other
> builtins aren't able to handle consistently we are looking for trouble.

But, then -funsafe-math-optimizations shouldn't be used.
Comment 9 Paolo Carlini 2013-07-25 16:15:14 UTC
Or maybe should be made a little weaker / safer? Are you 100% sure we are beating performancewise clang and icc on this?
Comment 10 Paolo Carlini 2013-07-25 17:27:25 UTC
Gaby, do you have an opinion on this? Irrespective of the long double issue, do you want me to re-enable (contra LWG 844) the pow(const complex<>&, int) overload in C++11 mode? If you think so, I can do it.
Comment 11 Gabriel Dos Reis 2013-07-25 18:02:56 UTC
(In reply to Paolo Carlini from comment #10)
> Gaby, do you have an opinion on this? Irrespective of the long double issue,
> do you want me to re-enable (contra LWG 844) the pow(const complex<>&, int)
> overload in C++11 mode? If you think so, I can do it.

I think pow(const complex<T>&,int) is lesser evil.  So, yes, I will support re-enabling it.  What would we run afoul, except the obvious standard letter?  Would we be computing something wrong?

-- Gaby
Comment 12 Paolo Carlini 2013-07-25 18:14:09 UTC
Agreed, let's do it.
Comment 13 Paolo Carlini 2013-07-25 21:24:41 UTC
Done.
Comment 14 Henner Sudek 2013-07-26 11:49:38 UTC
First i want to thank you all for the quick response and solving my problem. I just want to point out that std::pow(std::complex<long double>(0,0),1.) still returns nan. Maybe there is an way to unify the behavior of these functions?
Comment 15 Paolo Carlini 2013-07-26 11:53:29 UTC
However, there is *nothing* new about that. I still do believe there is a middle-end issue here, if only because clang and icc are fine.
Comment 16 Paolo Carlini 2013-07-26 12:20:50 UTC
Also, in practice, I think it's *very* hard to explain to the users why long double is so special, why the middle-end can't handle it in complete analogy with float and double. And since clang and icc are *already* doing it, apparently, you can't just tell them, vaguely, "it's very tough to implement" or "performance would be horrible" or something similar.
Comment 17 Marc Glisse 2013-07-26 12:58:01 UTC
(In reply to Henner Sudek from comment #14)
> First i want to thank you all for the quick response and solving my problem.
> I just want to point out that std::pow(std::complex<long double>(0,0),1.)
> still returns nan. Maybe there is an way to unify the behavior of these
> functions?

One thing that may be missing is a middle-end optimization that detects that the second argument to cpow is a simple constant (we have that for the real pow). But that wouldn't help for a runtime value of 1.

(In reply to Paolo Carlini from comment #16)
> Also, in practice, I think it's *very* hard to explain to the users why long
> double is so special, why the middle-end can't handle it in complete analogy
> with float and double. And since clang and icc are *already* doing it,
> apparently, you can't just tell them, vaguely, "it's very tough to
> implement" or "performance would be horrible" or something similar.

Did you benchmark both versions? That would help move the conversation...

I really don't think we want to give users too high expectations about -ffast-math. When you use complex numbers and the pow function, you know that along the negative real axis be dragons. Expecting -ffast-math to give standard results there is wrong, even if you may be lucky sometimes. We do a number of other transformations with fast-math that can be very wrong. Like we could have rounded the first argument of pow to a very small negative number before the call ;-)
Comment 18 Paolo Carlini 2013-07-26 13:23:14 UTC
Couple of clarifications: this doesn't go through cpow at all, the second argument isn't complex; this isn't -ffast-math, and in general in my experience whatever you throw at clang and icc (in fact, clang++ accepts the very same -f gcc accepts) like -Ofast, it never happens that exp(-inf) becomes nan (now I don't have the time to actually run benchmarks but since we are talking about compile-time constants I don't think we are beating them performance-wise on this); why long double is so special for exp(-inf), I don't think we can make a case for that, in terms of consistency or x87 crazyness (we are comparing also to icc not to a random compiler from people not really knowing the architecture).

But sorry I'm not going to contribute much more to this discussion, at the moment I don't think I'm able to make the compiler concretely better and it only would make me nervous.
Comment 19 Marc Glisse 2013-07-26 13:52:29 UTC
Ah, I think I understand now. For the exp function, getting an infinite argument violates the assumptions of fast-math, so anything goes (though of course it would be good to have the right result if it wasn't slower). But the user used the pow function, which even if it has a branch cut along the negative axis remains close to 0 around 0. The C++ standard does define pow in terms of exp and log, but I understand why it is surprising in this case.

Ok, forget I posted anything, I am not taking any side...

Ah, your reply arrived in the meantime:

(In reply to Paolo Carlini from comment #18)
> Couple of clarifications: this doesn't go through cpow at all, the second
> argument isn't complex;

Might have made sense to convert the second argument to a complex, call __builtin_cpow and let the middle-end deal with it.

> this isn't -ffast-math,

-funsafe-math-optimizations is pretty close to -ffast-math...

> and in general in my
> experience whatever you throw at clang and icc (in fact, clang++ accepts the
> very same -f gcc accepts) like -Ofast, it never happens that exp(-inf)
> becomes nan (now I don't have the time to actually run benchmarks

If someone does, that would help.

> but since
> we are talking about compile-time constants I don't think we are beating
> them performance-wise on this);

Don't we handle constant propagation using libmpc? That looks like a bug to me, we should indeed compute a constant in the middle-end.

> why long double is so special for exp(-inf),

It isn't long double that is special, it is float/double that get lucky ;-)
Though I admit that this surprising behavior is only ok if it comes with a significant speed gain.
Comment 20 Marc Glisse 2013-07-26 14:52:52 UTC
For reference, on my system, the performance ratio between a call to the libc expl and the inlined unsafe code is roughly 1.6.