[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Octave-bug-tracker] [bug #54302] power operator: 0.^0 results in NaN in
From: |
Dan Sebald |
Subject: |
[Octave-bug-tracker] [bug #54302] power operator: 0.^0 results in NaN in some case |
Date: |
Sat, 14 Jul 2018 04:55:09 -0400 (EDT) |
User-agent: |
Mozilla/5.0 (X11; Ubuntu; Linux x86_64; rv:55.0) Gecko/20100101 Firefox/55.0 |
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/