bugGNU Octave - Bugs: bug #54302, power operator: 0.^0 results in...

 
 

bug #54302: power operator: 0.^0 results in NaN in some case

Submitter:  None
Submitted:  Fri 13 Jul 2018 07:22:49 AM UTC
   
 
Category:  Octave Function Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Incorrect Result
Status:  Confirmed Assigned to:  None
Originator Name:  Originator Email:  -email is unavailable-
Open/Closed:  * Open Release:  * 4.4.0
Operating System:  * GNU/Linux Fixed Release:  None
Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Jump to the original submission

Fri 28 Feb 2020 06:20:47 PM UTC, comment #7: 

I know this issue report is old, but I just came across the issue again.  I wrote a long e-mail to the Maintainer's list with some of the some conclusions as were reached here.

Rik <rik5>
Group administrator
Mon 16 Jul 2018 03:44:34 AM UTC, comment #6: 

I was just thinking back, regarding the


 xpow (const Complex& a, double b)
 {
-  Complex result;
-
-  if (xisint (b))
-    result = std::pow (a, static_cast<int> (b));


that I removed in the patch, is it possibly the case that at one time these routines used


    result = pow (a, static_cast<int> (b));


and that was supposed to find the version of pow() I listed in the previous post, originating from oct-inttypes.cc?  That integer version might be slightly more accurate than std::pow(), but possibly a little slower...not drastically slow because it does use a fairly fast reduction by squaring the running product and bit-shifting b.  I would imagine that std::pow() is using logarithms or something similar to turn power into multiplication.

Dan Sebald <sebald>
Mon 16 Jul 2018 12:09:03 AM UTC, comment #5: 

I have made a modification that seems solve the problem.  I basically re-implemented all seven cases of std::pow and pow from the C++ library.  But that didn't fix the scalar operator results.  I then found a whole collection of scalar x scalar/matrix routines in xpow.cc and then replaced all of those with the new overloaded routine (zzpow ()).  Consequently, the casting to int cases I discarded; they seemed out-dated in the sense neither of the pow() functions have an integer exponent version.  There is


     double pow (Type1 base      , Type2 exponent);        // additional overloads


but the documentation says that most of those end up casting to a double...so it is like casting to an int just to be cast back to a double.  Maybe I didn't implement that correctly.  Anyway, it was producing errors so I took out the static_cast<int> which simplified some code.  See if you agree with what I changed there.

Here are a couple other important cases:


octave:1> diag([0 0]).^complex(0,0)
ans =

   NaN + NaNi   NaN + NaNi
   NaN + NaNi   NaN + NaNi

octave:2> zeros(2).^complex(0,0)
ans =

   NaN + NaNi   NaN + NaNi
   NaN + NaNi   NaN + NaNi


I'm a bit confused by all the combinations in xpow.cc versus the inline matrix implementations in mx-inlines.cc.  Do a


grep "pow (" */*/* -s


and one will see an awful lot of uses of std::pow().  So maybe those need a cursory review to confirm there aren't more x^0 issues, say for sparse and so on.

And then there is this from the oct-inttypes.cc file:


template <typename T>
octave_int<T>
pow (const octave_int<T>& a, const octave_int<T>& b)
{
  octave_int<T> retval;

  octave_int<T> zero = static_cast<T> (0);
  octave_int<T> one = static_cast<T> (1);

  if (b == zero || a == one)
    retval = one;
  else if (b < zero)
    {
      if (a == -one)
        retval = (b.value () % 2) ? a : one;
      else
        retval = zero;
    }
  else
    {
      octave_int<T> a_val = a;
      T b_val = b; // no need to do saturation on b

      retval = a;

      b_val -= 1;

      while (b_val != 0)
        {
          if (b_val & 1)
            retval = retval * a_val;

          b_val = b_val >> 1;

          if (b_val)
            a_val = a_val * a_val;
        }
    }

  return retval;
}


The above is a thorough consideration of all the base=1 and exponent=0 cases.  (Is it worth testing for base==1?  I mean, in most cases it will not be 1 and the test is just wasting a few cycles because the routine will produce a valid result if base=1.)  So, it's not like this issue hasn't been considered.  It's just that the std::pow routine has found extensive use throughout.  Then there is a separate powf () which may be repetitive and can be replaced by simple use of the overloaded pow () function.  Possibly some cruft here?

Lastly, I removed the following comment:


@@ -1096,19 +1104,6 @@ elem_xpow (const ComplexMatrix& a, const
 //
 //   * -> not needed.

-// FIXME: these functions need to be fixed so that things like
-//
-//   a = -1; b = [ 0, 0.5, 1 ]; r = a .^ b
-//
-// and
-//
-//   a = -1; b = [ 0, 0.5, 1 ]; for i = 1:3, r(i) = a .^ b(i), end
-//
-// produce identical results.  Also, it would be nice if -1^0.5
-// produced a pure imaginary result instead of a complex number with a
-// small real part.  But perhaps that's really a problem with the math
-// library...
-
 // -*- 1 -*-
 octave_value
 elem_xpow (double a, const NDArray& b)


I checked, and I found that two cases in the example given do in fact now match.  As for the (-realX)^I.5 comment, i.e., that this scenario should have the real portion forced to 0 (purely imaginary) I'm a bit hesitant to do that.  (We should be able to devise a short command to assess if the exponent has a fractional portion of exactly 0.5 if wanted, but the question is if we want that.)  My thinking is that the std::pow routine is using an algorithm to compute this value and if we force a certain subset of the values to be something slightly different than what the algorithm generates, it's as if we are creating a discontinuity, algorithm-wise.  That is, say, someone computes a series of numbers that converges to a -real value and looks at the results of std::pow(); there might be some subtle discontinue in some residual that catches the observer off guard.  To summarize, I'm thinking that if it is algorithmic in nature, probably best to leave it as is, but in the case of 0^0 it is an undefined input that we are simply assigning a more desired value to.  The user could certainly introduce logic as I described above to make (-realX)^I.5 purely imaginary.

(file #44562)

Dan Sebald <sebald>
Sat 14 Jul 2018 05:30:27 PM UTC, comment #4: 

The answer is yes:
Short answer is complex(0,0) is complex, 0+0j is not

octave:1>  0+0j
ans = 0
octave:2> whos
Variables in the current scope:

   Attr Name        Size                     Bytes  Class
   ==== ====        ====                     =====  =====
        ans         1x1                          8  double

Total is 1 element using 8 bytes

octave:3> complex(0,0)
ans =  0 + 0i
octave:4> whos
Variables in the current scope:

   Attr Name        Size                     Bytes  Class
   ==== ====        ====                     =====  =====
   c    ans         1x1                         16  double

Total is 1 element using 16 bytes

Michael Godfrey <godfrey>
Group Member
Sat 14 Jul 2018 08:55:09 AM UTC, comment #3: 

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.

Dan Sebald <sebald>
Sat 14 Jul 2018 06:12:52 AM UTC, comment #2: 

As a curious aside, is there a difference between the following?


octave:26> 0+0j
ans = 0
octave:27> complex(0,0)
ans =  0 + 0i


Dan Sebald <sebald>
Fri 13 Jul 2018 01:10:10 PM UTC, comment #1: 

I confirm that Matlab gives 1 in all the cases.

Marco Caliari <caliari>
Group Member
Fri 13 Jul 2018 07:22:49 AM UTC, original submission:  

Normally 0.^0 is 1, but in some case it results in NaN.


octave:495>  0 .^ 0
ans =  1
octave:496>  complex(0) .^ 0
ans =  1
octave:497>  0 .^ complex(0)
ans =  NaN + NaNi
octave:498>  complex(0) .^ complex(0)
ans =  NaN + NaNi



octave:502>  [1 2j;3j 0] .^ 0
ans =

   1   1
   1   1

octave:503>  [1 2j;3j 0] .^ [0 0;0 0]
ans =

   1   1
   1   1

octave:504>  [1 2j;3j 0] .^ [0 0]
ans =

     1 +   0i     1 +   0i
     1 +   0i   NaN - NaNi

octave:505>  [1 2j;3j 0] .^ [0;0]
ans =

     1 +   0i     1 +   0i
     1 +   0i   NaN - NaNi

octave:507>  [3j 0] .^ [0;0]
ans =

     1 +   0i   NaN - NaNi
     1 +   0i   NaN - NaNi

octave:508>  [2j;0] .^ [0 0]
ans =

     1 +   0i     1 +   0i
   NaN - NaNi   NaN - NaNi


Anonymous

 

(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attach Files:
   
   
Comment:
   

Attached Files

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by rik5 (Posted a comment)
  • -email is unavailable- added by godfrey (Posted a comment)
  • -email is unavailable- added by sebald (Posted a comment)
  • -email is unavailable- added by caliari (Posted a comment)
  • -email is unavailable- added by None (Submitted the item)
  •  

    There are 0 votes so far. Votes easily highlight which items people would like to see resolved in priority, independently of the priority of the item set by tracker managers.

    Only group members can vote.

     

    Follow 2 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2018-07-16 sebald Attached File- Added octave-zero_to_zero_power-djs2018jul15.patch, #44562
    2018-07-13 caliari StatusNone Confirmed

    Back to the top

    Powered by Savane 3.13-758e.
    Corresponding source code