[Bug libstdc++/13450] New: std::pow(std::complex<double>(-1.,0.),0.5) yields (NaN,0)

["pkienzle at nist dot gov" <gcc-bugzilla at gcc dot gnu dot org>]

In gcc/libstdc++-v3/std/std_complex.h between revision 1.5 and 1.6, the
definition of pow(cplx,real) was changed from:

  return exp(__y * log(__x));

to:

   if (__x.imag() == _Tp())
      return pow(__x.real(), __y);
 	 
   complex<_Tp> __t = log(__x);
   return polar(exp(__y * __t.real()), __y * __t.imag());

It should be testing for negative numbers as well:

   if (__x.imag() == _Tp() && __x.real() >= _Tp())
      ...


Similarly, pow(real,cplx) should also test for real<0.  

Instead of:

   return __x == _Tp()
      ? _Tp()
      : polar(pow(__x, __y.real()), __y.imag() * log(__x));

use something like:

   if (__x == _Tp()) 
      return _Tp();
   else if (__x > _Tp()) 
      return polar(pow(__x, __y.real()), __y.imag() * log(__x));
   else
      return pow(complex<_Tp>(__x,_Tp()), __y);


The following program demonstrates the problem:

#include <cmath>
#include <iostream>
#include <complex>

using namespace std;

typedef double _Tp;

complex<_Tp> mypow(complex<_Tp> __x, _Tp __y)
{
   if (__x.imag() == _Tp() && __x.real() >= _Tp())
      return pow(__x.real(), __y);
 	 
   complex<_Tp> __t = log(__x);
   return polar(exp(__y * __t.real()), __y * __t.imag());
}

complex<_Tp> mypow(_Tp __x, complex<_Tp> __y)
{
  if (__x == _Tp()) 
    return _Tp();
  else if (__x > _Tp())
    return polar(pow(__x, __y.real()), __y.imag() * log(__x));
  else
    return pow(complex<_Tp>(__x,_Tp()),__y);
}

typedef complex<double> cplx;
void test(double a, double b)
{
  cout << "a=" << a << ", b=" << b << endl;
  cout << "pow(-cplx,cplx)   =" << pow(cplx(a,0),cplx(b,0)) << endl;
  cout << "pow(-real,cplx)   =" << pow(a,cplx(b,0)) << endl;
  cout << "pow(-cplx,real)   =" << pow(cplx(a,0),b) << endl;
  cout << "mypow(-real,cplx) =" << mypow(a,cplx(b,0)) << endl;
  cout << "mypow(-cplx,real) =" << mypow(cplx(a,0),b) << endl;
}

int main(int argc, char *argv[])
{
  test(0,0.5);
  test(-1,0.5);
  test(-3.2,1.4);
  test(3.2,1.4);
  return 0;
}

-- 
           Summary: std::pow(std::complex<double>(-1.,0.),0.5) yields
                    (NaN,0)
           Product: gcc
           Version: 3.4.0
            Status: UNCONFIRMED
          Severity: normal
          Priority: P2
         Component: libstdc++
        AssignedTo: unassigned at gcc dot gnu dot org
        ReportedBy: pkienzle at nist dot gov
                CC: gcc-bugs at gcc dot gnu dot org


http://gcc.gnu.org/bugzilla/show_bug.cgi?id=13450