[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 #16 from Steve Kargl <sgk at troutmask dot apl.washington.edu> ---
On Mon, Apr 08, 2019 at 07:20:22PM +0000, redi at gcc dot gnu.org wrote:
> https://gcc.gnu.org/bugzilla/show_bug.cgi?id=89991
> 
> --- Comment #13 from Jonathan Wakely <redi at gcc dot gnu.org> ---
> (In reply to Steve Kargl from comment #10)
> > %  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)
> > 
> > This looks wrong.
> 
> I can't reproduce this, I get:
> 
>             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)
> 

My results are for i585-*-freebsd, which doesn't use glibc.
If Andrew is correct and a builtin is called, you might find
my results if you use -fno-builtins (check spelling).

Looking at ./libstdc++-v3/include/std/complex, one finds.

  // 26.2.8/13 sqrt(__z): Returns the complex square root of __z.
  //                     The branch cut is on the negative axis.
  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);
        }
      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);
        }
    }

Doesn't this gets the wrong answer for __y = -0 (as -0 < 0 is false)?