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

[Dan Sebald]

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

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

Additional Item Attachment:

File name: octave-zero_to_zero_power-djs2018jul15.patch Size:36 KB


    _______________________________________________________

Reply to this item at:

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

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