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)?