[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 #11 from Steve Kargl <sgk at troutmask dot apl.washington.edu> ---
On Mon, Apr 08, 2019 at 02:32:38PM +0000, 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.0000000000000000e+00)
>   cpow(z, 0.5) = CMPLX( 2.6283607659835831e-20, -4.2924388775825818e-04)
>       csqrt(z) = CMPLX( 0.0000000000000000e+00, -4.2924388775825818e-04)
> csqrt(conj(z)) = CMPLX( 0.0000000000000000e+00,  4.2924388775825818e-04)
> conj(csqrt(z)) = CMPLX( 0.0000000000000000e+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<class T> complex<T> sqrt(const complex<T>& 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.