From bug-octave-request at bevo dot che dot wisc dot edu Fri Dec 19 11:46:41 2003
Subject: Re: (-1).^.5 == NaN
From: Paul Kienzle <pkienzle at jazz dot ncnr dot nist dot gov>
To: "John W. Eaton" <jwe at bevo dot che dot wisc dot edu>
Cc: bug-octave at bevo dot che dot wisc dot edu
Date: Fri, 19 Dec 2003 12:46:11 -0500

> OK, how about the following change?  OTOH, I'm not sure that this
> should be "fixed" in Octave, but maybe it can eventually be changed
> back when the buggy libm function is no longer likely to be
> encountered.

How about fixing it in lo-mappers instead?  We can include a test
in configure for buggy pow since IRIX's libstdc++ is buggy, too.

...

> so if X is real, it calls pow(double,double).  Oops.  I think it
> should be smarter than that.  Can you please report this to the
> libc/libm maintainers as well, since it is apparently their bug?  Or
> at least it seems like a bug, but there may be a statement in some
> (stupid) standard somewhere that requires this behavior.

Reported.  The source is actually part of the gcc tree instead of
the glibc tree.

I did not report the bug to SGI yet.


Here is some code which demonstrates the problem:


#define USE_FIX

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

using namespace std;

typedef double _Tp;

#ifdef USE_FIX
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);
}
#endif

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;
#ifdef USE_FIX
  cout << "mypow(real,cplx) =" << mypow(a,cplx(b,0)) << endl;
  cout << "mypow(cplx,real) =" << mypow(cplx(a,0),b) << endl;
#endif
}

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;
}



-------------------------------------------------------------
Octave is freely available under the terms of the GNU GPL.

Octave's home on the web:  http://www.octave.org
How to fund new projects:  http://www.octave.org/funding.html
Subscription information:  http://www.octave.org/archive.html
-------------------------------------------------------------