element-wise exponent operator

[Rik]

jwe, maintainers,

I started trying to understand bug #57671 (Power operator operating on
complex numbers results in NaN when broadcasting) a week ago.  The problem
behavior is

([1;0]*j) .^ [1, 0]
ans =

   0.0000 + 1.0000i   1.0000 +      0i
        0 +      0i      NaN -    NaNi

But, if the matrices are full so that there is no broadcasting then Octave
works correctly

([1,1; 0,0]*j) .^ [1, 0; 1, 0]
ans =

   0 + 1i   1 + 0i
   0 + 0i   1 + 0i

After obtuse debugging sessions with gdb, I found that it is not
broadcasting per se which is the problem, but how Octave is calling
std::pow from the C++ library.  The bug manifests in many different ways. 
For example, besides broadcasting, real scalars behave differently than
complex scalars.

0.^0
ans = 1

But, for complex scalars,

complex (0,0) .^ complex (0,0)
ans =  NaN + NaNi

Bugs #57908, #57671, #57476, #54302, #52706 are all manifestations of the
same issue.

My first request is for an explanation of the convoluted mapping strategy
that Octave uses for operators.  There are tons of C++ files in
libinterp/operators/ which seem to map all of the various combinations of a
scalar, vector, or matrix object with another scalar, vector, or matrix
object.  And this structure is also echoed in liboctave/operators.  And
then there are exceptions to this structure such as
libinterp/corefcn/xpow.cc which slots in to the architecture as it gets
called from libinterp/operators, but isn't stored there.  Why aren't we
just mapping the '+' operator to a BLAS call?  Or maybe we are and I don't
understand how we get there.

For bug #57671, the code which ends up getting called is in xpow.cc case #10.

// -*- 10 -*-
octave_value
elem_xpow (const ComplexNDArray& a, const NDArray& b)
{
  dim_vector a_dims = a.dims ();
  dim_vector b_dims = b.dims ();

  if (a_dims != b_dims)
    {
      if (! is_valid_bsxfun ("operator .^", a_dims, b_dims))
        octave::err_nonconformant ("operator .^", a_dims, b_dims);

      return bsxfun_pow (a, b);
    }

  ComplexNDArray result (a_dims);

  for (octave_idx_type i = 0; i < a.numel (); i++)
    {
      octave_quit ();
      double btmp = b(i);
      if (xisint (btmp))
        result(i) = std::pow (a(i), static_cast<int> (btmp));
      else
        result(i) = std::pow (a(i), btmp);
    }

  return result;
}

For the non-broadcasting case, which works, the key is to note the check
for an integer power (xisint (btmp)).  If this is found it calls std::pow
from the C++ library with a different function signature.  When
broadcasting is used it calls bsxfun_pow and this ends up calling code from
liboctave.  In particular, the individual operation is in
liboctave/operators/mx-inlines.cc

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;

  for (size_t i = 0; i < n; i++)
    r[i] = pow (x[i], y[i]);
}

This calls std::pow with a function signature of two reals which results in
NaN + NaNi.

At a minimum, all the uses of std::pow in xpow.cc and mx-inlines.cc should
be checked to see if they need to install a test for an integer and then
perform the same cast as is done for case 10.

Handling Complex values will be a little more difficult.  The function
signatures available for pow are

template<class T> complex<T> pow (const complex<T>& x, int y);
template<class T> complex<T> pow (const complex<T>& x, const complex<T>& y);
template<class T> complex<T> pow (const complex<T>& x, const T& y);
template<class T> complex<T> pow (const T& x, const complex<T>& y);

There is no signature for pow (complex<real>, complex<int>) which means we
likely need to check first for whether the exponent can be narrowed to
real, and then checked subsequently for whether it is an integer.  As an
example,

[0, j] .^ complex (0,0)
ans =

   NaN + NaNi     1 +   0i

But, if the code in xpow.cc is modified as shown below

// -*- 11 -*-
octave_value
elem_xpow (const ComplexNDArray& a, const Complex& b)
{
  ComplexNDArray result (a.dims ());

  bool isint = (b.imag() == 0) && xisint (b.real ());
  int bint = static_cast<int> (b.real ());

  for (octave_idx_type i = 0; i < a.numel (); i++)
    {
      octave_quit ();
      if (isint)
        result(i) = std::pow (a(i), bint);
      else
        result(i) = std::pow (a(i), b);
    }

  return result;
}


[0, j] .^ complex (0,0)
ans =

   1   1

All of this will be tedious, because there are so many cases, which makes
one wonder whether the code in xpow.cc can't make better use of templates.

Regardless, once fixed we also need to add BIST tests for all the weird
corner cases.

--Rik