[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 #10 from Steve Kargl <sgk at troutmask dot apl.washington.edu> ---
On Mon, Apr 08, 2019 at 09:59:22AM +0000, redi at gcc dot gnu.org wrote:
> https://gcc.gnu.org/bugzilla/show_bug.cgi?id=89991
> 
> --- Comment #7 from Jonathan Wakely <redi at gcc dot gnu.org> ---
> I think it's allowed. The standards have very little to say about accuracy of
> any mathematical functions, and complex<double>(0, 0.0) == complex<double>(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 <complex.h>
#include <stdio.h>

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


mobile:kargl[210] cat a.cpp
#include <complex>
#include <iomanip>
#include <iostream>

int
main(int argc, char *argv[])
{
   std::complex<double> z, t0, t1, t2, t3;
   z  = std::complex<double>(-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.