[Bug libstdc++/50880] __complex_acosh() picks wrong complex branch

[kreckel at ginac dot de]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=50880

Richard B. Kreckel  changed:

   What|Removed |Added

   See Also||http://gcc.gnu.org/bugzilla
   ||/show_bug.cgi?id=50957

--- Comment #11 from Richard B. Kreckel  2011-11-02 
08:26:51 UTC ---
Paolo, I still intend to come forward with a patch for all these cases.
Unfortunately, I was distracted by what I've just filed as Bug 50957.

[Bug libstdc++/50880] __complex_acosh() picks wrong complex branch

[paolo.carlini at oracle dot com]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=50880

--- Comment #12 from Paolo Carlini  2011-11-02 
09:40:26 UTC ---
In my opinion BC2 is fine, I can take of applying it, if you still endorse it.

[Bug libstdc++/50880] __complex_acosh() picks wrong complex branch

[gdr at gcc dot gnu.org]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=50880

Gabriel Dos Reis  changed:

   What|Removed |Added

 CC||gdr at gcc dot gnu.org

--- Comment #13 from Gabriel Dos Reis  2011-11-02 
12:22:26 UTC ---
(In reply to comment #2)
> Thus, to understand and clarify why this has not been noticed so far, you are
> on a target which doesn't support in the underlying C library these complex
> functions, right? Because normally (eg, on Linux) these days we just forward 
> to
> __builtin_cacosh*, the code you are touching is just a "surrogate", a
> "fallback", which doesn't get right all the special cases, NaNs, infinity.
> 
> Anyway, a similar tweak would touch also the C++11 version in std::
> 
> Gaby, can you have a look to this, double check the patch? For your 
> convenience
> the surrounding code is:
> 
>   template
> std::complex<_Tp>
> __complex_acosh(const std::complex<_Tp>& __z)
> {
>   std::complex<_Tp> __t((__z.real() - __z.imag())
> * (__z.real() + __z.imag()) - _Tp(1.0),
> _Tp(2.0) * __z.real() * __z.imag());
>   __t = std::sqrt(__t);
> 
>   return std::log(__t + __z);
> }

As I observed elsewhere, the test should be on the sign, no comparison
against 0.0, so that signed zero is handled correctly.

[Bug libstdc++/50880] __complex_acosh() picks wrong complex branch

[gdr at gcc dot gnu.org]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=50880

--- Comment #14 from Gabriel Dos Reis  2011-11-02 
12:27:20 UTC ---
(In reply to comment #9)
> Created attachment 25654 [details]
> BC2

Since we are talking about branch cut and prespectiving
since zeros, I think we should avoid the 
the arithmetic z -/+ one, whee one is of a complex.
Rather the computation should be be directly on
the components.  This is to prevent signed zeros
to have their mutated.

[Bug libstdc++/50880] __complex_acosh() picks wrong complex branch

[paolo.carlini at oracle dot com]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=50880

--- Comment #15 from Paolo Carlini  2011-11-02 
12:31:09 UTC ---
Ok, thanks for your feedback Gaby. Indeed, I also wondered if we shouldn't work
with the components.

Richard, can you send a version of Kahan's algorithm rewritten in terms of real
and imag, I think we can all agree on that and resolve the PR for good.

[Bug libstdc++/50880] __complex_acosh() picks wrong complex branch

[paolo.carlini at oracle dot com]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=50880

--- Comment #16 from Paolo Carlini  2011-11-02 
12:44:06 UTC ---
Well, I guess this would be most of it:

  template
std::complex<_Tp>
__complex_acosh(const std::complex<_Tp>& __z)
{
  return _Tp(2.0) * std::log(std::sqrt(_Tp(0.5) * (__z + _Tp(1.0)))
 + std::sqrt(_Tp(0.5) * (__z - _Tp(1.0;
}

[Bug libstdc++/50880] __complex_acosh() picks wrong complex branch

[gdr at gcc dot gnu.org]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=50880

--- Comment #17 from Gabriel Dos Reis  2011-11-02 
12:48:23 UTC ---
(In reply to comment #16)
> Well, I guess this would be most of it:
> 
>   template
> std::complex<_Tp>
> __complex_acosh(const std::complex<_Tp>& __z)
> {
>   return _Tp(2.0) * std::log(std::sqrt(_Tp(0.5) * (__z + _Tp(1.0)))
>  + std::sqrt(_Tp(0.5) * (__z - _Tp(1.0;
> }

looks good -- hoping for log implementation to do the right thing.

[Bug libstdc++/50880] __complex_acosh() picks wrong complex branch

[gdr at gcc dot gnu.org]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=50880

--- Comment #18 from Gabriel Dos Reis  2011-11-02 
12:48:47 UTC ---
(In reply to comment #16)
> Well, I guess this would be most of it:
> 
>   template
> std::complex<_Tp>
> __complex_acosh(const std::complex<_Tp>& __z)
> {
>   return _Tp(2.0) * std::log(std::sqrt(_Tp(0.5) * (__z + _Tp(1.0)))
>  + std::sqrt(_Tp(0.5) * (__z - _Tp(1.0;
> }

looks good -- hoping for log implementation to do the right thing.

[Bug libstdc++/50880] __complex_acosh() picks wrong complex branch

[paolo at gcc dot gnu.org]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=50880

--- Comment #19 from paolo at gcc dot gnu.org  
2011-11-02 18:43:31 UTC ---
Author: paolo
Date: Wed Nov  2 18:43:26 2011
New Revision: 180787

URL: http://gcc.gnu.org/viewcvs?root=gcc&view=rev&rev=180787
Log:
2011-11-02  Richard B. Kreckel  
Paolo Carlini  

PR libstdc++/50880
* include/std/complex (__complex_acosh): Fix in a better way,
use Kahan's formula.
* include/tr1/complex (__complex_acosh): Likewise.

Modified:
trunk/libstdc++-v3/include/std/complex
trunk/libstdc++-v3/include/tr1/complex

[Bug libstdc++/50880] __complex_acosh() picks wrong complex branch

[paolo at gcc dot gnu.org]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=50880

--- Comment #20 from paolo at gcc dot gnu.org  
2011-11-02 18:43:46 UTC ---
Author: paolo
Date: Wed Nov  2 18:43:42 2011
New Revision: 180788

URL: http://gcc.gnu.org/viewcvs?root=gcc&view=rev&rev=180788
Log:
2011-11-02  Richard B. Kreckel  
Paolo Carlini  

PR libstdc++/50880
* include/std/complex (__complex_acosh): Fix in a better way,
use Kahan's formula.
* include/tr1/complex (__complex_acosh): Likewise.

Modified:
trunk/libstdc++-v3/ChangeLog

[Bug libstdc++/50880] __complex_acosh() picks wrong complex branch

[paolo at gcc dot gnu.org]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=50880

--- Comment #21 from paolo at gcc dot gnu.org  
2011-11-02 21:54:29 UTC ---
Author: paolo
Date: Wed Nov  2 21:54:24 2011
New Revision: 180804

URL: http://gcc.gnu.org/viewcvs?root=gcc&view=rev&rev=180804
Log:
2011-11-02  Richard B. Kreckel  
Paolo Carlini  

PR libstdc++/50880
* include/std/complex (__complex_acosh): Fix in a better way,
use Kahan's formula.
* include/tr1/complex (__complex_acosh): Likewise.

Added:
branches/gcc-4_6-branch/libstdc++-v3/testsuite/26_numerics/complex/50880.cc
   
branches/gcc-4_6-branch/libstdc++-v3/testsuite/tr1/8_c_compatibility/complex/50880.cc
Modified:
branches/gcc-4_6-branch/libstdc++-v3/ChangeLog
branches/gcc-4_6-branch/libstdc++-v3/include/std/complex
branches/gcc-4_6-branch/libstdc++-v3/include/tr1/complex

[Bug libstdc++/50880] __complex_acosh() picks wrong complex branch

[paolo.carlini at oracle dot com]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=50880

Paolo Carlini  changed:

   What|Removed |Added

 Status|NEW |RESOLVED
 Resolution||FIXED

--- Comment #22 from Paolo Carlini  2011-11-02 
21:56:15 UTC ---
Fixed for 4.6.3 and mainline.

[Bug libstdc++/50880] __complex_acosh() picks wrong complex branch

[kreckel at ginac dot de]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=50880

--- Comment #23 from Richard B. Kreckel  2011-11-03 
23:57:55 UTC ---
(In reply to comment #16)
> Well, I guess this would be most of it:
> 
>   template
> std::complex<_Tp>
> __complex_acosh(const std::complex<_Tp>& __z)
> {
>   return _Tp(2.0) * std::log(std::sqrt(_Tp(0.5) * (__z + _Tp(1.0)))
>  + std::sqrt(_Tp(0.5) * (__z - _Tp(1.0;
> }

[Sorry for my temporary absence.]

For future reference, some final remarks:
* Yes, that is a good implementation for this "fallback".
  It relies on __z - _Tp(1.0) not "mutating" the sign of __z's imag part.
* If the complex log does not do the right thing, all is lost anyways
  (besides, __complex_asinh relies on it, too).
* My patch BC1 was flawed. It contains code trying to work around a ctor doing
  multiplication (fixed in PR48760)
* My patch BC2 was flawed for the same reason: __z - __one doesn't preserve
  the sign of __z's imag part.

Looks good. Thanks!

[Bug libstdc++/50880] __complex_acosh() picks wrong complex branch

[paolo.carlini at oracle dot com]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=50880

--- Comment #24 from Paolo Carlini  2011-11-04 
00:51:49 UTC ---
Thanks Richard for double checking!

[Bug libstdc++/50880] __complex_acosh() picks wrong complex branch

[paolo.carlini at oracle dot com]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=50880

--- Comment #25 from Paolo Carlini  2011-11-04 
00:53:29 UTC ---
By the way, if isn't clear already, I would be *really* curious to know which
specific targets by now can't just enable the builtins, eg, their libc doesn't
provide the C99 complex functions.

[Bug libstdc++/50880] __complex_acosh() picks wrong complex branch

[kreckel at ginac dot de]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=50880

--- Comment #26 from Richard B. Kreckel  2011-11-04 
08:17:20 UTC ---
(In reply to comment #25)
> By the way, if isn't clear already, I would be *really* curious to know which
> specific targets by now can't just enable the builtins, eg, their libc doesn't
> provide the C99 complex functions.

Sadly, the libc functions aren't always correct either. But that is another
story: http://sourceware.org/bugzilla/show_bug.cgi?id=13305

[Bug libstdc++/50880] __complex_acosh() picks wrong complex branch

[kreckel at ginac dot de]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=50880

--- Comment #1 from Richard B. Kreckel  2011-10-27 
07:12:12 UTC ---
Created attachment 25623
  --> http://gcc.gnu.org/bugzilla/attachment.cgi?id=25623
patch to fix the bug

[Bug libstdc++/50880] __complex_acosh() picks wrong complex branch

[paolo.carlini at oracle dot com]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=50880

Paolo Carlini  changed:

   What|Removed |Added

   Keywords|wrong-code  |
 CC||g...@integrable-solutions.ne
   ||t

--- Comment #2 from Paolo Carlini  2011-10-27 
09:24:26 UTC ---
Thus, to understand and clarify why this has not been noticed so far, you are
on a target which doesn't support in the underlying C library these complex
functions, right? Because normally (eg, on Linux) these days we just forward to
__builtin_cacosh*, the code you are touching is just a "surrogate", a
"fallback", which doesn't get right all the special cases, NaNs, infinity.

Anyway, a similar tweak would touch also the C++11 version in std::

Gaby, can you have a look to this, double check the patch? For your convenience
the surrounding code is:

  template
std::complex<_Tp>
__complex_acosh(const std::complex<_Tp>& __z)
{
  std::complex<_Tp> __t((__z.real() - __z.imag())
* (__z.real() + __z.imag()) - _Tp(1.0),
_Tp(2.0) * __z.real() * __z.imag());
  __t = std::sqrt(__t);

  return std::log(__t + __z);
}

[Bug libstdc++/50880] __complex_acosh() picks wrong complex branch

[paolo.carlini at oracle dot com]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=50880

Paolo Carlini  changed:

   What|Removed |Added

 Status|UNCONFIRMED |NEW
   Last reconfirmed||2011-10-27
   Target Milestone|--- |4.6.3
 Ever Confirmed|0   |1

--- Comment #3 from Paolo Carlini  2011-10-27 
09:49:16 UTC ---
Ok, checked. Will commit the fix momentarily (4.6.3 too when the branch
reopens)

[Bug libstdc++/50880] __complex_acosh() picks wrong complex branch

[paolo at gcc dot gnu.org]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=50880

--- Comment #4 from paolo at gcc dot gnu.org  
2011-10-27 11:00:30 UTC ---
Author: paolo
Date: Thu Oct 27 11:00:25 2011
New Revision: 180563

URL: http://gcc.gnu.org/viewcvs?root=gcc&view=rev&rev=180563
Log:
2011-10-27  Richard B. Kreckel  
Paolo Carlini  

PR libstdc++/50880
* include/std/complex (__complex_acosh): Fix for __z.real() < 0.
* include/tr1/complex (__complex_acosh): Likewise.
* testsuite/26_numerics/complex/50880.cc: New.
* testsuite/tr1/8_c_compatibility/complex/50880.cc: Likewise.

Added:
trunk/libstdc++-v3/testsuite/26_numerics/complex/50880.cc
trunk/libstdc++-v3/testsuite/tr1/8_c_compatibility/complex/50880.cc
Modified:
trunk/libstdc++-v3/ChangeLog
trunk/libstdc++-v3/include/std/complex
trunk/libstdc++-v3/include/tr1/complex

[Bug libstdc++/50880] __complex_acosh() picks wrong complex branch

[kreckel at ginac dot de]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=50880

--- Comment #5 from Richard B. Kreckel  2011-10-28 
07:06:57 UTC ---
On 10/27/2011 11:24 AM, paolo.carlini at oracle dot com wrote:
> Thus, to understand and clarify why this has not been noticed so far, you are
> on a target which doesn't support in the underlying C library these complex
> functions, right? Because normally (eg, on Linux) these days we just forward 
> to
> __builtin_cacosh*, the code you are touching is just a "surrogate", a
> "fallback", which doesn't get right all the special cases, NaNs, infinity.

Well, I didn't "notice" it. Searching for bugs involving branch cut 
positions of inverse trig functions in a range of FOSS projects is a pet 
project of mine.  ;-)

BTW: I noticed that my patch tests if (__z.real() < -0.0), which is 
correct, but the sign of zero looks eccentric and might potentially 
confuse future readers. I suggest removing it. It doesn't matter at all, 
anyhow.

[Bug libstdc++/50880] __complex_acosh() picks wrong complex branch

[paolo.carlini at oracle dot com]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=50880

--- Comment #6 from Paolo Carlini  2011-10-28 
09:10:09 UTC ---
Indeed you are right about the sign, in terms at least of consistency with the
rest of the fallback implementations which already have got quite a number of
comparisons with zero with no special attention to its signedness (like '<
_Tp()' or '> _Tp()'). I had already noticed that. As soon as I find a bit of
time, we can also *consistently over all those cases* use __builtin_signbit, as
suggested by Gaby elsewhere. I have to double check with the middle-end people
that it doesn't generate library calls for the patch to be neat.

[Bug libstdc++/50880] __complex_acosh() picks wrong complex branch

[kreckel at ginac dot de]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=50880

--- Comment #7 from Richard B. Kreckel  2011-10-28 
21:52:08 UTC ---
> As soon as I find a bit of
> time, we can also *consistently over all those cases* use __builtin_signbit, 
> as
> suggested by Gaby elsewhere. I have to double check with the middle-end people
> that it doesn't generate library calls for the patch to be neat.

We also better double check whether the results stay correct.

Thinking of it... Big Ooops!

It turns out the patch makes it much better but still not entirely correct. On
systems that feature a signed zero, it gives incorrect results when either
* __z.real() is -0.0 and signbit(__z.imag())
* __z.real() < -1.0 and __z.imag() is -0.0

The first problem can be fixed by using signbit instead of -0.0, as you
proposed, but the second needs more correction. The patch BC1 I'm attaching
should fix these cases, too.

But this is starting to look cumbersome. We are trying to construct a formula
that is continuous except on the branch cut defined in C99. Why not just use
the formula mentioned in CLTL2 [0] like in patch BC2 that I'm attaching? (The
branch cuts of inverse trig functinos in C99 and Common Lisp are the same.)

[0]

[Bug libstdc++/50880] __complex_acosh() picks wrong complex branch

[kreckel at ginac dot de]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=50880

--- Comment #8 from Richard B. Kreckel  2011-10-28 
21:53:30 UTC ---
Created attachment 25653
  --> http://gcc.gnu.org/bugzilla/attachment.cgi?id=25653
BC1

[Bug libstdc++/50880] __complex_acosh() picks wrong complex branch

[kreckel at ginac dot de]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=50880

--- Comment #9 from Richard B. Kreckel  2011-10-28 
21:54:07 UTC ---
Created attachment 25654
  --> http://gcc.gnu.org/bugzilla/attachment.cgi?id=25654
BC2

[Bug libstdc++/50880] __complex_acosh() picks wrong complex branch

[paolo.carlini at oracle dot com]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=50880

--- Comment #10 from Paolo Carlini  2011-10-28 
22:09:57 UTC ---
Richard, I have no problems with BC2. This is code I wrote rather quickly a few
years ago, adapting it from glibc, essentially, and then each year that went
by, fewer and fewer systems used and tested it, because underlying C99 libc
support is becoming often available (eg, for sure Linux and Darwin). Thus,
please sleep on this, let's wait for comments from other people, and say, in a
week or so we'll finalize the code for 4.7. Thanks again!