[Octave-bug-tracker] [bug #54302] power operator: 0.^0 results in NaN in some case

[Dan Sebald]

Follow-up Comment #3, bug #54302 (project octave):

The ^ operator and power() function call end up using the same op_el_pow
construction.  It's a long series of defines that builds a table of operators
for which lookup is done to get the function to call.

I can see in the file mx-inlines.cc that at least in the case of matrices and
vectors (may not be the same for scalar operations, don't know), the
individual power computations are done with the routines:


mx_inline_pow (size_t n, R *r, const X *x, const Y *y)
mx_inline_pow (size_t n, R *r, const X *x, Y y)
mx_inline_pow (size_t n, R *r, X x, const Y *y)


These inline routines use the std::pow() function which is a complex function.
 [There is a comment within that suggests the compiler is left to choose
pow(), but the "using" directive might force std::pow().  The other option
would be cmath pow(double, double).  There is also a FIXME questioning this. 
Should clear this up.]  The reference for std::pow is here (both C98 and
C++11):

http://www.cplusplus.com/reference/complex/pow/

and the reference for cmath pow is here:

http://www.cplusplus.com/reference/cmath/pow/

In both cases, the base of 0 and exponent of 0 (either double or complex)
results in a returned value which is application specific, whereas the error
flags that are thrown are more well-defined.  So, it seems to me either the
options are to check for the exponent being 0 prior, or check for error flags
after calling pow() and then decipher.  Probably the former.

This creates a problem, though, in the sense that these inline routines are
defined as templates and we can't generally do


if (y[i] == 0)
  r[i] = 1;


The following fails as well for cases where Y and R are not complex, e.g.,
just "float":


diff --git a/liboctave/operators/mx-inlines.cc
b/liboctave/operators/mx-inlines.cc
--- a/liboctave/operators/mx-inlines.cc
+++ b/liboctave/operators/mx-inlines.cc
@@ -404,18 +405,19 @@ DEFMINMAXSPEC (double, mx_inline_xmax, >
 DEFMINMAXSPEC (float, mx_inline_xmin, <=)
 DEFMINMAXSPEC (float, mx_inline_xmax, >=)
 
-// FIXME: Is this comment correct anymore?  It seems like std::pow is
chosen.
-// Let the compiler decide which pow to use, whichever best matches the
-// arguments provided.
-
 template <typename R, typename X, typename Y>
 inline void
 mx_inline_pow (size_t n, R *r, const X *x, const Y *y)
 {
   using std::pow;
+  std::complex<double> zero(0, 0);
+  static std::complex<double> one(1, 0);
 
   for (size_t i = 0; i < n; i++)
-    r[i] = pow (x[i], y[i]);
+    if (y[i] == zero)
+      r[i] = one;
+    else
+      r[i] = pow (x[i], y[i]);
 }
 
 template <typename R, typename X, typename Y>
@@ -423,9 +425,15 @@ inline void
 mx_inline_pow (size_t n, R *r, const X *x, Y y)
 {
   using std::pow;
+  static std::complex<double> zero(0, 0);
+  static std::complex<double> one(1, 0);
 
-  for (size_t i = 0; i < n; i++)
-    r[i] = pow (x[i], y);
+  if (y == zero)
+    for (size_t i = 0; i < n; i++)
+      r[i] = one;
+  else
+    for (size_t i = 0; i < n; i++)
+      r[i] = pow (x[i], y);
 }
 
 template <typename R, typename X, typename Y>
@@ -433,9 +441,14 @@ inline void
 mx_inline_pow (size_t n, R *r, X x, const Y *y)
 {
   using std::pow;
+  std::complex<double> zero(0, 0);
+  static std::complex<double> one(1, 0);
 
   for (size_t i = 0; i < n; i++)
-    r[i] = pow (x, y[i]);
+    if (y[i] == zero)
+      r[i] = one;
+    else
+      r[i] = pow (x, y[i]);
 }
 
 // Arbitrary function appliers.


So, there is a not-so-easy task of either expanding out all the different
scenarios based on R, X and Y typenames.  Or, write the two-operand
mx_inline_XYZ routines to have another two inputs, call them ZERO and ONE as
follows


mx_inline_pow (size_t n, R *r, const X *x, Y y, Y ZERO, R ONE)


where the y is compared against ZERO and if equal then r = ONE.  But that
means all two-argument operator functions have to change, even though the
inlining would completely ignore those last two values if they don't appear
within the operator function.

This isn't a fun change.

    _______________________________________________________

Reply to this item at:

  <http://savannah.gnu.org/bugs/?54302>

_______________________________________________
  Message sent via Savannah
  https://savannah.gnu.org/