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

2019-04-09 Thread sgk at troutmask dot apl.washington.edu
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

2019-04-09 Thread redi at gcc dot gnu.org
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

2019-04-08 Thread kargl at gcc dot gnu.org
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

2019-04-08 Thread sgk at troutmask dot apl.washington.edu
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

2019-04-08 Thread kargl at gcc dot gnu.org
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

2019-04-08 Thread redi at gcc dot gnu.org
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

2019-04-08 Thread 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 +, 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

2019-04-08 Thread sgk at troutmask dot apl.washington.edu
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

2019-04-08 Thread pinskia at gcc dot gnu.org
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

2019-04-08 Thread sgk at troutmask dot apl.washington.edu
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

2019-04-08 Thread pinskia at gcc dot gnu.org
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

2019-04-08 Thread redi at gcc dot gnu.org
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

2019-04-08 Thread redi at gcc dot gnu.org
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

2019-04-08 Thread redi at gcc dot gnu.org
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

2019-04-08 Thread sgk at troutmask dot apl.washington.edu
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

2019-04-08 Thread sgk at troutmask dot apl.washington.edu
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

2019-04-08 Thread redi at gcc dot gnu.org
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

2019-04-08 Thread pinskia at gcc dot gnu.org
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

2019-04-08 Thread redi at gcc dot gnu.org
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

2019-04-08 Thread pinskia at gcc dot gnu.org
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

2019-04-08 Thread rguenth at gcc dot gnu.org
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

2019-04-08 Thread marxin at gcc dot gnu.org
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

2019-04-06 Thread t.sprodowski at web dot de
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

2019-04-05 Thread kargl at gcc dot gnu.org
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

2019-04-05 Thread t.sprodowski at web dot de
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.