[Bug libstdc++/89991] Complex numbers: Calculation of imaginary part is not correct

[sgk at troutmask dot apl.washington.edu]

https://gcc.gnu.org/bugzilla/show_bug.cgi?id=89991

--- Comment #18 from Steve Kargl <sgk at troutmask dot apl.washington.edu> ---
On Mon, Apr 08, 2019 at 08:03:36PM +0000, pinskia at gcc dot gnu.org wrote:
> https://gcc.gnu.org/bugzilla/show_bug.cgi?id=89991
> 
> --- Comment #17 from Andrew Pinski <pinskia at gcc dot gnu.org> ---
> >Doesn't this gets the wrong answer for __y = -0 (as -0 < 0 is false)?
> 
> No, you missed this part:
> //                     The branch cut is on the negative axis.

No, I didn't miss that part.

> So maybe the bug is inside FreeBSD and Window's libm.  Glibc fixed the branch
> cuts issues back in 2012 for csqrt but the other OS's did not change theirs.

For the C++ code in comment, on x86_64-*-freebsd.

% g++8 -o z a.cpp -lm && ./z
            z = (-1.84250315177824e-07,-0)
   pow(z,0.5) = (2.62836076598358e-20,-0.000429243887758258)
      sqrt(z) = (0,0.000429243887758258)
sqrt(conj(z)) = (0,0.000429243887758258)
conj(sqrt(z)) = (0,-0.000429243887758258)

The last two lines are definitely wrong.

troutmask:sgk[209] nm z | grep csqrt
troutmask:sgk[210] nm z | grep sqrt
000000000040156b W _ZSt14__complex_sqrtIdESt7complexIT_ERKS2_
000000000040143d W _ZSt4sqrtIdESt7complexIT_ERKS2_
                 U sqrt@@FBSD_1.0

There is no reference to csqrt in the exectuable.  If I change
/usr/local/lib/gcc8/include/c++/complex to use copysign
to account for __y = -0 like

  template<typename _Tp>
    complex<_Tp>
    __complex_sqrt(const complex<_Tp>& __z)
    {
      _Tp __x = __z.real();
      _Tp __y = __z.imag();

      if (__x == _Tp())
        {
          _Tp __t = sqrt(abs(__y) / 2);
//          return complex<_Tp>(__t, __y < _Tp() ? -__t : __t);
          return complex<_Tp>(__t, copysign(__t, __y));
        }
      else
        {
          _Tp __t = sqrt(2 * (std::abs(__z) + abs(__x)));
          _Tp __u = __t / 2;
//          return __x > _Tp()
//            ? complex<_Tp>(__u, __y / __t)
//            : complex<_Tp>(abs(__y) / __t, __y < _Tp() ? -__u : __u);
          return __x > _Tp()
            ? complex<_Tp>(__u, __y / __t)
            : complex<_Tp>(abs(__y) / __t, copysign(__u, __y));
        }
    }


The C++ code in comment #10 gives

 g++8 -o z a.cpp -lm && ./z
            z = (-1.84250315177824e-07,-0)
   pow(z,0.5) = (2.62836076598358e-20,-0.000429243887758258)
      sqrt(z) = (0,-0.000429243887758258)
sqrt(conj(z)) = (0,0.000429243887758258)
conj(sqrt(z)) = (0,0.000429243887758258)

The correct answer.  QED.