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

[kargl at gcc dot gnu.org]

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

--- Comment #19 from kargl at gcc dot gnu.org ---
(In reply to Steve Kargl from comment #18)
> 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.

BTW, if change /usr/local/lib/gcc8/include/c++/complex back to
no using copysign(), and instead change

#if _GLIBCXX_USE_C99_COMPLEX
  inline __complex__ float
  __complex_sqrt(__complex__ float __z) { return __builtin_csqrtf(__z); }

to 

#if _GLIBCXX_USE_C99_COMPLEX || SOMETHING_UGLY
  inline __complex__ float
  __complex_sqrt(__complex__ float __z) { return __builtin_csqrtf(__z); }

and do 

 g++8 -DSOMETHING_UGLY -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)

I get the expected.  So, if you're on a system that has
_GLIBCXX_USE_C99_COMPLEX, you won't see the bug.

It is likely that everywhere that a construct of the
form __y < _Tp() ? -__u : __u appear, it needs to use
copysign.