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

https://gcc.gnu.org/bugzilla/show_bug.cgi?id=89991 --- Comment #25 from Steve Kargl --- On Tue, Apr 09, 2019 at 08:24:29PM +, redi at gcc dot gnu.org wrote: > https://gcc.gnu.org/bugzilla/show_bug.cgi?id=89991 > > --- Comment #24 from Jonathan Wakely --- > Thanks for the patch, I'll test it fully tomorrow. > I think the patch for complex_sqrt() is correct. The one for complex_pow(), I think accidently works for OP, but is likely broken for some general regions of the complex plane. > I'll open a separate bug for the FreeBSD issue. We could use more fine-grained > configure checks so that most C99 math functions are enabled, even if some of > the complex ones are missing. libgfortran has c99_functions.c that implements missing C99 math functions when configure cannot find one. The implementations are likely to be fairly direct without much optimization, worrying about exceptional casea, or even tested extensively.

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

https://gcc.gnu.org/bugzilla/show_bug.cgi?id=89991 --- Comment #24 from Jonathan Wakely --- Thanks for the patch, I'll test it fully tomorrow. I'll open a separate bug for the FreeBSD issue. We could use more fine-grained configure checks so that most C99 math functions are enabled, even if some of the complex ones are missing.

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

https://gcc.gnu.org/bugzilla/show_bug.cgi?id=89991 --- Comment #23 from kargl at gcc dot gnu.org --- Created attachment 46105 --> https://gcc.gnu.org/bugzilla/attachment.cgi?id=46105=edit fix g++ problem with pow(z,0.5) where imag(z) = -0. This patch has only been tested with the original test provided by the reporter.

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

https://gcc.gnu.org/bugzilla/show_bug.cgi?id=89991 --- Comment #22 from Steve Kargl --- On Mon, Apr 08, 2019 at 10:17:00PM +, redi at gcc dot gnu.org wrote: > (In reply to kargl from comment #19) > > I get the expected. So, if you're on a system that has > > _GLIBCXX_USE_C99_COMPLEX, you won't see the bug. > > Wow, why isn't libstdc++ using the C99 functions on FreeBSD? > Because it is all or nothing. See comment #8 and #9 in https://gcc.gnu.org/bugzilla/show_bug.cgi?id=89125 If you're too busy to look at PR89125 the upshot is that FreeBSD is missing ccoshl, ccosl, cexpl, csinhl, csinl, ctanhl, and ctanl. I have BSD licensed versions of ccoshl, ccosl, cexpl, csinhl, and csinl, but testing on FreeBSD ran into what I consider to be a very bad bug in clang (FreeBSD system compiler).

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

https://gcc.gnu.org/bugzilla/show_bug.cgi?id=89991 --- Comment #21 from kargl at gcc dot gnu.org --- Created attachment 46102 --> https://gcc.gnu.org/bugzilla/attachment.cgi?id=46102=edit fix g++ problem with sqrt(z) where z is complex and imag(z) = -0

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

https://gcc.gnu.org/bugzilla/show_bug.cgi?id=89991 --- Comment #20 from Jonathan Wakely --- (In reply to Steve Kargl from comment #16) > If Andrew is correct and a builtin is called, Wasn't that me, not Andrew? > you might find > my results if you use -fno-builtins (check spelling). No, same results. Calling __builtin_csqrt doesn't necessarily mean GCC evaluates it. Without optimisation it still generates a call to csqrt from libc. > 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 > 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)? Yes, but that code shouldn't be used for modern targets ... (In reply to kargl from comment #19) > I get the expected. So, if you're on a system that has > _GLIBCXX_USE_C99_COMPLEX, you won't see the bug. Wow, why isn't libstdc++ using the C99 functions on FreeBSD? I'll have to look into that. > It is likely that everywhere that a construct of the > form __y < _Tp() ? -__u : __u appear, it needs to use > copysign. That won't always work, because the generic functions should really only get used when _Tp is not one of float, double or long double. And in that case there might be no copysign for the type. For float, double and long double we should be using the libc routines. So the bug is that FreeBSD isn't using them. The _original_ bug report is for std::pow though, and on Ubuntu, which does use glibc. Comment 9 needs more analysis.

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

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 +, pinskia at gcc dot gnu.org wrote: > > https://gcc.gnu.org/bugzilla/show_bug.cgi?id=89991 > > > > --- Comment #17 from Andrew Pinski --- > > >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 > 0040156b W _ZSt14__complex_sqrtIdESt7complexIT_ERKS2_ > 0040143d 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 > 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.

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

https://gcc.gnu.org/bugzilla/show_bug.cgi?id=89991 --- Comment #18 from Steve Kargl --- On Mon, Apr 08, 2019 at 08:03:36PM +, pinskia at gcc dot gnu.org wrote: > https://gcc.gnu.org/bugzilla/show_bug.cgi?id=89991 > > --- Comment #17 from Andrew Pinski --- > >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 0040156b W _ZSt14__complex_sqrtIdESt7complexIT_ERKS2_ 0040143d 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 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.

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

https://gcc.gnu.org/bugzilla/show_bug.cgi?id=89991 --- Comment #17 from Andrew Pinski --- >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. 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.

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

https://gcc.gnu.org/bugzilla/show_bug.cgi?id=89991 --- Comment #16 from Steve Kargl --- On Mon, Apr 08, 2019 at 07:20:22PM +, redi at gcc dot gnu.org wrote: > https://gcc.gnu.org/bugzilla/show_bug.cgi?id=89991 > > --- Comment #13 from Jonathan Wakely --- > (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 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)?

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

https://gcc.gnu.org/bugzilla/show_bug.cgi?id=89991 --- Comment #15 from Andrew Pinski --- (In reply to Jonathan Wakely from comment #13) > (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 bet now comes to the fact there have been improvements to glibc which changed the behavior here Also I used the wrong term, it is the branch cut that is the issue. Most of the branch cuts were fixed in glibc in 2012; though there might have been some fixed later on.

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

https://gcc.gnu.org/bugzilla/show_bug.cgi?id=89991 --- Comment #14 from Jonathan Wakely --- Which is unsurprising because std::sqrt(z) just calls __builtin_csqrt(z.__rep())

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

https://gcc.gnu.org/bugzilla/show_bug.cgi?id=89991 --- Comment #13 from Jonathan Wakely --- (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)

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

https://gcc.gnu.org/bugzilla/show_bug.cgi?id=89991 --- Comment #12 from Jonathan Wakely --- (In reply to Steve Kargl from comment #11) > unless [Note: ...] is non-normative text. That's exactly what it is. But we can still aim to meet the intended behaviour.

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

https://gcc.gnu.org/bugzilla/show_bug.cgi?id=89991 --- Comment #11 from Steve Kargl --- On Mon, Apr 08, 2019 at 02:32:38PM +, sgk at troutmask dot apl.washington.edu wrote: > > I don't have a copy of the C++ standard, so take this specualtion. > pow(z,0.5) is equivalent to sqrt(z). From the C standard, one has > > conj(csqrt(z)) = csqrt(conj(z)). > > g++ does not enforce this when the imaginary part is -0; > while gcc does. (code snipped) > % gcc8 -o z c.c -lm && ./z > z = CMPLX(-1.8425031517782417e-07, -0.e+00) > cpow(z, 0.5) = CMPLX( 2.6283607659835831e-20, -4.2924388775825818e-04) > csqrt(z) = CMPLX( 0.e+00, -4.2924388775825818e-04) > csqrt(conj(z)) = CMPLX( 0.e+00, 4.2924388775825818e-04) > conj(csqrt(z)) = CMPLX( 0.e+00, 4.2924388775825818e-04) (code snipped) > % 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. It is wrong. From n4810.pdf, page 1102, template complex sqrt(const complex& x); Returns: The complex square root of x, in the range of the right half-plane. [Note: The semantics of this function are intended to be the same in C++ as they are for csqrt in C. -- end note] unless [Note: ...] is non-normative text.

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

https://gcc.gnu.org/bugzilla/show_bug.cgi?id=89991 --- Comment #10 from Steve Kargl --- On Mon, Apr 08, 2019 at 09:59:22AM +, redi at gcc dot gnu.org wrote: > https://gcc.gnu.org/bugzilla/show_bug.cgi?id=89991 > > --- Comment #7 from Jonathan Wakely --- > I think it's allowed. The standards have very little to say about accuracy of > any mathematical functions, and complex(0, 0.0) == complex(0, > -0.0) is true according to the standard, because +0.0 == -0.0 is true. > I don't have a copy of the C++ standard, so take this specualtion. pow(z,0.5) is equivalent to sqrt(z). From the C standard, one has conj(csqrt(z)) = csqrt(conj(z)). g++ does not enforce this when the imaginary part is -0; while gcc does. % cat c.c #include #include int main(void) { double complex z, t0, t1, t2, t3; z = CMPLX(-1.8425031517782417e-07, -0.0); t0 = cpow(z, 0.5); t1 = csqrt(z); t2 = csqrt(conj(z)); t3 = conj(csqrt(z)); printf(" z = CMPLX(% .16le, % .16le)\n", creal(z), cimag(z)); printf(" cpow(z, 0.5) = CMPLX(% .16le, % .16le)\n", creal(t0), cimag(t0)); printf(" csqrt(z) = CMPLX(% .16le, % .16le)\n", creal(t1), cimag(t1)); printf("csqrt(conj(z)) = CMPLX(% .16le, % .16le)\n", creal(t2), cimag(t2)); printf("conj(csqrt(z)) = CMPLX(% .16le, % .16le)\n", creal(t3), cimag(t3)); return 0; } % gcc8 -o z c.c -lm && ./z z = CMPLX(-1.8425031517782417e-07, -0.e+00) cpow(z, 0.5) = CMPLX( 2.6283607659835831e-20, -4.2924388775825818e-04) csqrt(z) = CMPLX( 0.e+00, -4.2924388775825818e-04) csqrt(conj(z)) = CMPLX( 0.e+00, 4.2924388775825818e-04) conj(csqrt(z)) = CMPLX( 0.e+00, 4.2924388775825818e-04) mobile:kargl[210] cat a.cpp #include #include #include int main(int argc, char *argv[]) { std::complex z, t0, t1, t2, t3; z = std::complex(-1.8425031517782417e-07, -0.0); t0 = std::pow(z, 0.5); t1 = std::sqrt(z); t2 = std::sqrt(std::conj(z)); t3 = std::conj(std::sqrt(z)); std::cout << "z = " << std::setprecision(15) << z << std::endl; std::cout << " pow(z,0.5) = " << std::setprecision(15) << t0 << std::endl; std::cout << " sqrt(z) = " << std::setprecision(15) << t1 << std::endl; std::cout << "sqrt(conj(z)) = " << std::setprecision(15) << t2 << std::endl; std::cout << "conj(sqrt(z)) = " << std::setprecision(15) << t3 << std::endl; return 0; } % 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.

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

https://gcc.gnu.org/bugzilla/show_bug.cgi?id=89991 --- Comment #9 from Jonathan Wakely --- (In reply to Richard Biener from comment #5) > The issue is > > std::pow (__x=..., __y=@0x7fffdcb8: 0.5) > at /home/space/rguenther/install/gcc-9.0/include/c++/9.0.1/complex:1027 > (gdb) l > 1022{ > 1023#if ! _GLIBCXX_USE_C99_COMPLEX > 1024 if (__x == _Tp()) > 1025return _Tp(); > 1026#endif > 1027 if (__x.imag() == _Tp() && __x.real() > _Tp()) > 1028return pow(__x.real(), __y); > > where __x.imag () == _Tp() says true for -0.0 == 0.0. This means > std::pow will return the same values for r + -0.0i and r + 0.0i, > not sure if that is allowed by the C++ standard. But __x.real() > _Tp() is false here, so that branch isn't taken anyway. Instead the pow(val, 0.5) result comes from: _Complex double val = -1.8425031517782417e-07 + -0.0 * I; _Complex double t = __builtin_clog(val); double rho = exp(0.5 * __real__ t); double theta = 0.5 * __imag__ t; _Complex result = rho * cos(theta) + rho * sin(theta) * I; __builtin_printf("%f\n", __imag__ result);

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

https://gcc.gnu.org/bugzilla/show_bug.cgi?id=89991 --- Comment #8 from Andrew Pinski --- Also isn't it true that this is just a different quadrant of the solution? That is the answer is correct but which quadrant being selected is different? That is (a^0.5) actually has two answers where the imaginary part can be positive or negative? That is they are conjugate of each other.

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

https://gcc.gnu.org/bugzilla/show_bug.cgi?id=89991 --- Comment #7 from Jonathan Wakely --- I think it's allowed. The standards have very little to say about accuracy of any mathematical functions, and complex(0, 0.0) == complex(0, -0.0) is true according to the standard, because +0.0 == -0.0 is true.

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

https://gcc.gnu.org/bugzilla/show_bug.cgi?id=89991 --- Comment #6 from Andrew Pinski --- (In reply to Richard Biener from comment #5) > The issue is > > std::pow (__x=..., __y=@0x7fffdcb8: 0.5) > at /home/space/rguenther/install/gcc-9.0/include/c++/9.0.1/complex:1027 > (gdb) l > 1022{ > 1023#if ! _GLIBCXX_USE_C99_COMPLEX > 1024 if (__x == _Tp()) > 1025return _Tp(); > 1026#endif > 1027 if (__x.imag() == _Tp() && __x.real() > _Tp()) > 1028return pow(__x.real(), __y); > > where __x.imag () == _Tp() says true for -0.0 == 0.0. This means > std::pow will return the same values for r + -0.0i and r + 0.0i, > not sure if that is allowed by the C++ standard. If it does not allow it, then adding copysign is needed.

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

https://gcc.gnu.org/bugzilla/show_bug.cgi?id=89991 Richard Biener changed: What|Removed |Added Keywords||wrong-code --- Comment #5 from Richard Biener --- The issue is std::pow (__x=..., __y=@0x7fffdcb8: 0.5) at /home/space/rguenther/install/gcc-9.0/include/c++/9.0.1/complex:1027 (gdb) l 1022{ 1023#if ! _GLIBCXX_USE_C99_COMPLEX 1024 if (__x == _Tp()) 1025return _Tp(); 1026#endif 1027 if (__x.imag() == _Tp() && __x.real() > _Tp()) 1028return pow(__x.real(), __y); where __x.imag () == _Tp() says true for -0.0 == 0.0. This means std::pow will return the same values for r + -0.0i and r + 0.0i, not sure if that is allowed by the C++ standard.

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

https://gcc.gnu.org/bugzilla/show_bug.cgi?id=89991 Martin Liška changed: What|Removed |Added Status|UNCONFIRMED |NEW Last reconfirmed||2019-04-08 CC||marxin at gcc dot gnu.org Ever confirmed|0 |1 --- Comment #4 from Martin Liška --- I see the expected result when replacing '-0.0' with '0.0'. Well, negative zero should be equal to the positive one according to standard.

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

https://gcc.gnu.org/bugzilla/show_bug.cgi?id=89991 --- Comment #3 from t.sprodowski at web dot de --- Octave 4.2.2: ans = 2.6284e-20 + 4.2924e-04i

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

https://gcc.gnu.org/bugzilla/show_bug.cgi?id=89991 kargl at gcc dot gnu.org changed: What|Removed |Added CC||kargl at gcc dot gnu.org --- Comment #2 from kargl at gcc dot gnu.org --- (In reply to t.sprodowski from comment #0) > Following calculation of the complex number leads to a wrong imaginary part: > > > #include > #include > #include > > int main(int argc, char *argv[]) > { > std::complex val = std::complex(-1.8425031517782417e-07, > -0.0); > std::complex testExp = std::pow(val, 0.5); > std::cout << "textExp: " << std::setprecision(30) << testExp << std::endl; > return 0; > } > > Result is: > (2.6283607659835830609796003783e-20,-0.000429243887758258178214548772544), > but it should be > (2.628360765983583e-20, 0.0004292438877582582), obtained from Visual Studio, > MATLAB and Octave. > What version of Octave. I get >> z = complex(-1.8425031517782417e-07, -0.0) z = -0.0018425 - 0.000i >> z**0.5 ans = 2.6284e-20 - 4.2924e-04i which agrees with clang++ version 7.0.1 (and apparently g++ which I haven't tested).

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

https://gcc.gnu.org/bugzilla/show_bug.cgi?id=89991 --- Comment #1 from t.sprodowski at web dot de --- Created attachment 46095 --> https://gcc.gnu.org/bugzilla/attachment.cgi?id=46095=edit Source file Source file to illustrate this bug.