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

On 19-Dec-2003, Paul Kienzle <pkienzle at jazz dot ncnr dot nist dot gov> wrote:

| Debian octave-2.1.52
| 
|    octave:9>  (-1).^.5
|    ans = NaN
| 
| IRIX octave-2.1.40
| 
|    octave:9>  (-1).^.5
|    ans = NaN
| 
| This is arguably a libm problem rather than an Octave
| problem, as the following code demonstrates:
| 
| #include <cmath>
| #include <iostream>
| #include <complex>
| 
| int main(int argc, char *argv[])
| {
|   std::cout << "real^real=" << std::pow(-1.,0.5) << std::endl;
|   std::cout << "real^cplx=" << std::pow(-1.,std::complex<double>(0.5,0)) << std::endl;
|   std::cout << "cplx^real=" << std::pow(std::complex<double>(-1.,0),0.5) << std::endl;
|   std::cout << "cplx^cplx=" << std::pow(std::complex<double>(-1.,0),std::complex<double>(0.5,0)) << std::endl;
|   return 0;
| }
| 
| Debian libm:
| 
|    real^real=nan
|    real^cplx=(nan,nan)
|    cplx^real=(nan,0)
|    cplx^cplx=(6.12303e-17,1)
| 
| IRIX libm:
| 
|   real^real=nan
|   real^cplx=(nan,nan)
|   cplx^real=(6.12323e-17,1)
|   cplx^cplx=(6.12323e-17,1)

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.

The <complex> header file for g++ 3.3 on Debian contains

  template<typename _Tp>
    complex<_Tp>
    pow(const complex<_Tp>& __x, const _Tp& __y)
    {
      if (__x.imag() == _Tp())
        return pow(__x.real(), __y);

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

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.

Thanks,

jwe


2003-12-19  John W. Eaton  <jwe at bevo dot che dot wisc dot edu>

	* xpow.cc (xpow (double, double)): Avoid apparent GNU libm bug.


Index: src/xpow.cc
===================================================================
RCS file: /usr/local/cvsroot/octave/src/xpow.cc,v
retrieving revision 1.35
diff -u -r1.35 xpow.cc
--- src/xpow.cc	26 Nov 2003 07:02:42 -0000	1.35
+++ src/xpow.cc	19 Dec 2003 15:12:51 -0000
 at @ -70,8 +70,13 @@
 {
   if (a < 0.0 && static_cast<int> (b) != b)
     {
+      // XXX FIXME XXX -- avoid apparent GNU libm bug by converting
+      // A and B to complex instead of just A.
+
       Complex atmp (a);
-      return pow (atmp, b);
+      Complex btmp (b);
+
+      return pow (atmp, btmp);
     }
   else
     return pow (a, b);



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