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

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

[kreckel at ginac dot de]

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

             Bug #: 50880
           Summary: __complex_acosh() picks wrong complex branch
    Classification: Unclassified
           Product: gcc
           Version: unknown
            Status: UNCONFIRMED
          Severity: normal
          Priority: P3
         Component: libstdc++
        AssignedTo: unassigned@gcc.gnu.org
        ReportedBy: kreckel@ginac.de


According to the C++ standard, acosh(complex<T>) should behave just like C99's
cacosh(T complex) function. There, the branch cut is "at values less than 1
along the real axis" and the "range of a half-strip of non-negative values
along the real axis and in the interval [-i*pi,+i*pi] along the imaginary
axis."

The implementation in tr1/complex gets this wrong. The result returned by
__complex_acosh() are all wrong in the lower complex plain. It can be easily
fixed. I'm attaching a patch.

[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 <kreckel at ginac dot de> 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 <paolo.carlini at oracle dot com> changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
           Keywords|wrong-code                  |
                 CC|                            |gdr@integrable-solutions.ne
                   |                            |t

--- Comment #2 from Paolo Carlini <paolo.carlini at oracle dot com> 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<typename _Tp>
    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 <paolo.carlini at oracle dot com> 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 <paolo.carlini at oracle dot com> 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 <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  <kreckel@ginac.de>
        Paolo Carlini  <paolo.carlini@oracle.com>

    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 <kreckel at ginac dot de> 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 <paolo.carlini at oracle dot com> 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 <kreckel at ginac dot de> 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]
<http://www-prod-gif.supelec.fr/docs/cltl/clm/node129.html#SECTION001653000000000000000>

[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 <kreckel at ginac dot de> 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

[kreckel at ginac dot de]

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

--- Comment #8 from Richard B. Kreckel <kreckel at ginac dot de> 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

[paolo.carlini at oracle dot com]

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

--- Comment #10 from Paolo Carlini <paolo.carlini at oracle dot com> 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!

[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 <kreckel at ginac dot de> changed:

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

--- Comment #11 from Richard B. Kreckel <kreckel at ginac dot de> 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 <paolo.carlini at oracle dot com> 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 <gdr at gcc dot gnu.org> changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
                 CC|                            |gdr at gcc dot gnu.org

--- Comment #13 from Gabriel Dos Reis <gdr at gcc dot gnu.org> 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<typename _Tp>
>     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 <gdr at gcc dot gnu.org> 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<T>.
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 <paolo.carlini at oracle dot com> 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 <paolo.carlini at oracle dot com> 2011-11-02 12:44:06 UTC ---
Well, I guess this would be most of it:

  template<typename _Tp>
    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 <gdr at gcc dot gnu.org> 2011-11-02 12:48:23 UTC ---
(In reply to comment #16)
> Well, I guess this would be most of it:
> 
>   template<typename _Tp>
>     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 <gdr at gcc dot gnu.org> 2011-11-02 12:48:47 UTC ---
(In reply to comment #16)
> Well, I guess this would be most of it:
> 
>   template<typename _Tp>
>     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 <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  <kreckel@ginac.de>
        Paolo Carlini  <paolo.carlini@oracle.com>

    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 <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  <kreckel@ginac.de>
        Paolo Carlini  <paolo.carlini@oracle.com>

    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 <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  <kreckel@ginac.de>
        Paolo Carlini  <paolo.carlini@oracle.com>

    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 <paolo.carlini at oracle dot com> changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
             Status|NEW                         |RESOLVED
         Resolution|                            |FIXED

--- Comment #22 from Paolo Carlini <paolo.carlini at oracle dot com> 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 <kreckel at ginac dot de> 2011-11-03 23:57:55 UTC ---
(In reply to comment #16)
> Well, I guess this would be most of it:
> 
>   template<typename _Tp>
>     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 <paolo.carlini at oracle dot com> 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 <paolo.carlini at oracle dot com> 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 <kreckel at ginac dot de> 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