exporting patch:
# HG changeset patch
# User Daniel J Sebald <daniel.sebald@ieee.org>
# Date 1531698800 18000
#      Sun Jul 15 18:53:20 2018 -0500
# Node ID 479788924c68b69fa05a61cf4d6c5514a9afede7
# Parent  6fb384c855c1f97b17f921ca4d7d9a23adbd1fd6
Make undefined 0^0 be 1 for all the combinations of variable types (bug #54302)

* xpow.cc: Replace all std::pow() by zzpow() which will set 0^0 to 1.  Remove
  all cases where a float or double which is cast to an int is fed into zzpow()
  and simply use the original float value.  The two cases where the exponent
  is -1 and normal division is used remain.  Remove FIXME comments about two
  methods that should produce same result (they did already) and raising
  negative real values to I.5 power being purely imaginary.
  (xpow (const ComplexDiagMatrix& a, const Complex& b)): In this case, if the
  exponent is 0, return an identity DiagMatrix rather than a ComplexDiagMatrix,
  making this consistent with xpow (const ComplexMatrix& a, double b).

* mx-inlines.cc: Drop the FIXME about std::pow being chosen.
  (zzpow): Replicate all the possible std::pow() and cmath.h pow() overloads
  that can exist in the C++11 compiler but make them do check for an exponent
  of 0, thereby makeing 0^0 be 1 rather than application specific.
  (mx_inline_pow): In all overload cases, drop "using std::pow" and replace
  pow() with zzpow().

diff --git a/libinterp/corefcn/xpow.cc b/libinterp/corefcn/xpow.cc
--- a/libinterp/corefcn/xpow.cc
+++ b/libinterp/corefcn/xpow.cc
@@ -100,10 +100,10 @@ xpow (double a, double b)
     {
       Complex atmp (a);
 
-      return std::pow (atmp, b);
+      return zzpow (atmp, b);
     }
   else
-    retval = std::pow (a, b);
+    retval = zzpow (a, b);
 
   return retval;
 }
@@ -131,9 +131,9 @@ xpow (double a, const Matrix& b)
         {
           Complex elt = lambda(i);
           if (std::imag (elt) == 0.0)
-            lambda(i) = std::pow (a, std::real (elt));
+            lambda(i) = zzpow (a, std::real (elt));
           else
-            lambda(i) = std::pow (a, elt);
+            lambda(i) = zzpow (a, elt);
         }
       ComplexDiagMatrix D (lambda);
 
@@ -155,7 +155,7 @@ xpow (double a, const Matrix& b)
 octave_value
 xpow (double a, const Complex& b)
 {
-  Complex result = std::pow (a, b);
+  Complex result = zzpow (a, b);
   return result;
 }
 
@@ -182,9 +182,9 @@ xpow (double a, const ComplexMatrix& b)
         {
           Complex elt = lambda(i);
           if (std::imag (elt) == 0.0)
-            lambda(i) = std::pow (a, std::real (elt));
+            lambda(i) = zzpow (a, std::real (elt));
           else
-            lambda(i) = std::pow (a, elt);
+            lambda(i) = zzpow (a, elt);
         }
       ComplexDiagMatrix D (lambda);
 
@@ -269,7 +269,7 @@ xpow (const Matrix& a, double b)
           ComplexMatrix Q (a_eig.right_eigenvectors ());
 
           for (octave_idx_type i = 0; i < nr; i++)
-            lambda(i) = std::pow (lambda(i), b);
+            lambda(i) = zzpow (lambda(i), b);
 
           ComplexDiagMatrix D (lambda);
 
@@ -300,14 +300,14 @@ xpow (const DiagMatrix& a, double b)
     {
       DiagMatrix r (nr, nc);
       for (octave_idx_type i = 0; i < nc; i++)
-        r.dgelem (i) = std::pow (a.dgelem (i), b);
+        r.dgelem (i) = zzpow (a.dgelem (i), b);
       retval = r;
     }
   else
     {
       ComplexDiagMatrix r (nr, nc);
       for (octave_idx_type i = 0; i < nc; i++)
-        r.dgelem (i) = std::pow (static_cast<Complex> (a.dgelem (i)), b);
+        r.dgelem (i) = zzpow (static_cast<Complex> (a.dgelem (i)), b);
       retval = r;
     }
 
@@ -344,7 +344,7 @@ xpow (const Matrix& a, const Complex& b)
       ComplexMatrix Q (a_eig.right_eigenvectors ());
 
       for (octave_idx_type i = 0; i < nr; i++)
-        lambda(i) = std::pow (lambda(i), b);
+        lambda(i) = zzpow (lambda(i), b);
 
       ComplexDiagMatrix D (lambda);
 
@@ -362,13 +362,7 @@ xpow (const Matrix& a, const Complex& b)
 octave_value
 xpow (const Complex& a, double b)
 {
-  Complex result;
-
-  if (xisint (b))
-    result = std::pow (a, static_cast<int> (b));
-  else
-    result = std::pow (a, b);
-
+  Complex result = zzpow (a, b);
   return result;
 }
 
@@ -395,9 +389,9 @@ xpow (const Complex& a, const Matrix& b)
         {
           Complex elt = lambda(i);
           if (std::imag (elt) == 0.0)
-            lambda(i) = std::pow (a, std::real (elt));
+            lambda(i) = zzpow (a, std::real (elt));
           else
-            lambda(i) = std::pow (a, elt);
+            lambda(i) = zzpow (a, elt);
         }
       ComplexDiagMatrix D (lambda);
 
@@ -415,8 +409,7 @@ xpow (const Complex& a, const Matrix& b)
 octave_value
 xpow (const Complex& a, const Complex& b)
 {
-  Complex result;
-  result = std::pow (a, b);
+  Complex result = zzpow (a, b);
   return result;
 }
 
@@ -443,9 +436,9 @@ xpow (const Complex& a, const ComplexMat
         {
           Complex elt = lambda(i);
           if (std::imag (elt) == 0.0)
-            lambda(i) = std::pow (a, std::real (elt));
+            lambda(i) = zzpow (a, std::real (elt));
           else
-            lambda(i) = std::pow (a, elt);
+            lambda(i) = zzpow (a, elt);
         }
       ComplexDiagMatrix D (lambda);
 
@@ -530,7 +523,7 @@ xpow (const ComplexMatrix& a, double b)
           ComplexMatrix Q (a_eig.right_eigenvectors ());
 
           for (octave_idx_type i = 0; i < nr; i++)
-            lambda(i) = std::pow (lambda(i), b);
+            lambda(i) = zzpow (lambda(i), b);
 
           ComplexDiagMatrix D (lambda);
 
@@ -565,7 +558,7 @@ xpow (const ComplexMatrix& a, const Comp
       ComplexMatrix Q (a_eig.right_eigenvectors ());
 
       for (octave_idx_type i = 0; i < nr; i++)
-        lambda(i) = std::pow (lambda(i), b);
+        lambda(i) = zzpow (lambda(i), b);
 
       ComplexDiagMatrix D (lambda);
 
@@ -591,10 +584,15 @@ xpow (const ComplexDiagMatrix& a, const 
   if (nr == 0 || nc == 0 || nr != nc)
     err_nonsquare_matrix ();
 
-  ComplexDiagMatrix r (nr, nc);
-  for (octave_idx_type i = 0; i < nc; i++)
-    r(i, i) = std::pow (a(i, i), b);
-  retval = r;
+  if (b == Complex (0, 0))
+    retval = DiagMatrix (nr, nr, 1.0);
+  else
+    {
+      ComplexDiagMatrix r (nr, nc);
+      for (octave_idx_type i = 0; i < nc; i++)
+        r(i, i) = zzpow (a(i, i), b);
+      retval = r;
+    }
 
   return retval;
 }
@@ -627,19 +625,6 @@ xpow (const DiagMatrix& a, const Complex
 //
 //   * -> 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 Matrix& b)
@@ -660,7 +645,7 @@ elem_xpow (double a, const Matrix& b)
         for (octave_idx_type i = 0; i < nr; i++)
           {
             octave_quit ();
-            result(i, j) = std::pow (atmp, b(i, j));
+            result(i, j) = zzpow (atmp, b(i, j));
           }
 
       retval = result;
@@ -673,7 +658,7 @@ elem_xpow (double a, const Matrix& b)
         for (octave_idx_type i = 0; i < nr; i++)
           {
             octave_quit ();
-            result(i, j) = std::pow (a, b(i, j));
+            result(i, j) = zzpow (a, b(i, j));
           }
 
       retval = result;
@@ -696,7 +681,7 @@ elem_xpow (double a, const ComplexMatrix
     for (octave_idx_type i = 0; i < nr; i++)
       {
         octave_quit ();
-        result(i, j) = std::pow (atmp, b(i, j));
+        result(i, j) = zzpow (atmp, b(i, j));
       }
 
   return result;
@@ -722,8 +707,8 @@ elem_xpow (double a, const Range& r)
       Matrix result (1, n);
       if (same_sign (r.base (), r.inc ()))
         {
-          double base = std::pow (a, r.base ());
-          double inc = std::pow (a, r.inc ());
+          double base = zzpow (a, r.base ());
+          double inc = zzpow (a, r.inc ());
           result(0) = base;
           for (octave_idx_type i = 1; i < n; i++)
             result(i) = (base *= inc);
@@ -731,8 +716,8 @@ elem_xpow (double a, const Range& r)
       else
         {
           // Don't use Range::limit () here.
-          double limit = std::pow (a, r.base () + (n-1) * r.inc ());
-          double inc = std::pow (a, -r.inc ());
+          double limit = zzpow (a, r.base () + (n-1) * r.inc ());
+          double inc = zzpow (a, -r.inc ());
           result(n-1) = limit;
           for (octave_idx_type i = n-2; i >= 0; i--)
             result(i) = (limit *= inc);
@@ -766,7 +751,7 @@ elem_xpow (const Matrix& a, double b)
 
             Complex atmp (a(i, j));
 
-            result(i, j) = std::pow (atmp, b);
+            result(i, j) = zzpow (atmp, b);
           }
 
       retval = result;
@@ -779,7 +764,7 @@ elem_xpow (const Matrix& a, double b)
         for (octave_idx_type i = 0; i < nr; i++)
           {
             octave_quit ();
-            result(i, j) = std::pow (a(i, j), b);
+            result(i, j) = zzpow (a(i, j), b);
           }
 
       retval = result;
@@ -829,7 +814,7 @@ done:
             octave_quit ();
             Complex atmp (a(i, j));
             Complex btmp (b(i, j));
-            complex_result(i, j) = std::pow (atmp, btmp);
+            complex_result(i, j) = zzpow (atmp, btmp);
           }
 
       retval = complex_result;
@@ -842,7 +827,7 @@ done:
         for (octave_idx_type i = 0; i < nr; i++)
           {
             octave_quit ();
-            result(i, j) = std::pow (a(i, j), b(i, j));
+            result(i, j) = zzpow (a(i, j), b(i, j));
           }
 
       retval = result;
@@ -864,7 +849,7 @@ elem_xpow (const Matrix& a, const Comple
     for (octave_idx_type i = 0; i < nr; i++)
       {
         octave_quit ();
-        result(i, j) = std::pow (Complex (a(i, j)), b);
+        result(i, j) = zzpow (Complex (a(i, j)), b);
       }
 
   return result;
@@ -889,7 +874,7 @@ elem_xpow (const Matrix& a, const Comple
     for (octave_idx_type i = 0; i < nr; i++)
       {
         octave_quit ();
-        result(i, j) = std::pow (Complex (a(i, j)), b(i, j));
+        result(i, j) = zzpow (Complex (a(i, j)), b(i, j));
       }
 
   return result;
@@ -908,11 +893,7 @@ elem_xpow (const Complex& a, const Matri
     for (octave_idx_type i = 0; i < nr; i++)
       {
         octave_quit ();
-        double btmp = b(i, j);
-        if (xisint (btmp))
-          result(i, j) = std::pow (a, static_cast<int> (btmp));
-        else
-          result(i, j) = std::pow (a, btmp);
+        result(i, j) = zzpow (a, b(i, j));
       }
 
   return result;
@@ -931,7 +912,7 @@ elem_xpow (const Complex& a, const Compl
     for (octave_idx_type i = 0; i < nr; i++)
       {
         octave_quit ();
-        result(i, j) = std::pow (a, b(i, j));
+        result(i, j) = zzpow (a, b(i, j));
       }
 
   return result;
@@ -952,8 +933,8 @@ elem_xpow (const Complex& a, const Range
 
       if (same_sign (r.base (), r.inc ()))
         {
-          Complex base = std::pow (a, r.base ());
-          Complex inc = std::pow (a, r.inc ());
+          Complex base = zzpow (a, r.base ());
+          Complex inc = zzpow (a, r.inc ());
           result(0) = base;
           for (octave_idx_type i = 1; i < n; i++)
             result(i) = (base *= inc);
@@ -961,8 +942,8 @@ elem_xpow (const Complex& a, const Range
       else
         {
           // Don't use Range::limit () here.
-          Complex limit = std::pow (a, r.base () + (n-1) * r.inc ());
-          Complex inc = std::pow (a, -r.inc ());
+          Complex limit = zzpow (a, r.base () + (n-1) * r.inc ());
+          Complex inc = zzpow (a, -r.inc ());
           result(n-1) = limit;
           for (octave_idx_type i = n-2; i >= 0; i--)
             result(i) = (limit *= inc);
@@ -985,25 +966,12 @@ elem_xpow (const ComplexMatrix& a, doubl
 
   ComplexMatrix result (nr, nc);
 
-  if (xisint (b))
-    {
-      int bint = static_cast<int> (b);
-      for (octave_idx_type j = 0; j < nc; j++)
-        for (octave_idx_type i = 0; i < nr; i++)
-          {
-            octave_quit ();
-            result(i, j) = std::pow (a(i, j), bint);
-          }
-    }
-  else
-    {
-      for (octave_idx_type j = 0; j < nc; j++)
-        for (octave_idx_type i = 0; i < nr; i++)
-          {
-            octave_quit ();
-            result(i, j) = std::pow (a(i, j), b);
-          }
-    }
+  for (octave_idx_type j = 0; j < nc; j++)
+    for (octave_idx_type i = 0; i < nr; i++)
+      {
+        octave_quit ();
+        result(i, j) = zzpow (a(i, j), b);
+      }
 
   return result;
 }
@@ -1027,11 +995,7 @@ elem_xpow (const ComplexMatrix& a, const
     for (octave_idx_type i = 0; i < nr; i++)
       {
         octave_quit ();
-        double btmp = b(i, j);
-        if (xisint (btmp))
-          result(i, j) = std::pow (a(i, j), static_cast<int> (btmp));
-        else
-          result(i, j) = std::pow (a(i, j), btmp);
+        result(i, j) = zzpow (a(i, j), b(i, j));
       }
 
   return result;
@@ -1050,7 +1014,7 @@ elem_xpow (const ComplexMatrix& a, const
     for (octave_idx_type i = 0; i < nr; i++)
       {
         octave_quit ();
-        result(i, j) = std::pow (a(i, j), b);
+        result(i, j) = zzpow (a(i, j), b);
       }
 
   return result;
@@ -1075,7 +1039,7 @@ elem_xpow (const ComplexMatrix& a, const
     for (octave_idx_type i = 0; i < nr; i++)
       {
         octave_quit ();
-        result(i, j) = std::pow (a(i, j), b(i, j));
+        result(i, j) = zzpow (a(i, j), b(i, j));
       }
 
   return result;
@@ -1096,19 +1060,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)
@@ -1122,7 +1073,7 @@ elem_xpow (double a, const NDArray& b)
       for (octave_idx_type i = 0; i < b.numel (); i++)
         {
           octave_quit ();
-          result(i) = std::pow (atmp, b(i));
+          result(i) = zzpow (atmp, b(i));
         }
 
       retval = result;
@@ -1133,7 +1084,7 @@ elem_xpow (double a, const NDArray& b)
       for (octave_idx_type i = 0; i < b.numel (); i++)
         {
           octave_quit ();
-          result(i) = std::pow (a, b(i));
+          result(i) = zzpow (a, b(i));
         }
 
       retval = result;
@@ -1151,7 +1102,7 @@ elem_xpow (double a, const ComplexNDArra
   for (octave_idx_type i = 0; i < b.numel (); i++)
     {
       octave_quit ();
-      result(i) = std::pow (a, b(i));
+      result(i) = zzpow (a, b(i));
     }
 
   return result;
@@ -1173,7 +1124,7 @@ elem_xpow (const NDArray& a, double b)
             {
               octave_quit ();
               Complex atmp (a(i));
-              result(i) = std::pow (atmp, b);
+              result(i) = zzpow (atmp, b);
             }
 
           retval = result;
@@ -1184,7 +1135,7 @@ elem_xpow (const NDArray& a, double b)
           for (octave_idx_type i = 0; i < a.numel (); i++)
             {
               octave_quit ();
-              result(i) = std::pow (a(i), b);
+              result(i) = zzpow (a(i), b);
             }
 
           retval = result;
@@ -1215,7 +1166,7 @@ elem_xpow (const NDArray& a, double b)
           for (octave_idx_type i = 0; i < a.numel (); i++)
             {
               octave_quit ();
-              result(i) = std::pow (a(i), ib);
+              result(i) = zzpow (a(i), ib);
             }
         }
 
@@ -1274,7 +1225,7 @@ done:
         {
           octave_quit ();
           Complex atmp (a(i));
-          complex_result(i) = std::pow (atmp, b(i));
+          complex_result(i) = zzpow (atmp, b(i));
         }
 
       retval = complex_result;
@@ -1286,7 +1237,7 @@ done:
       for (octave_idx_type i = 0; i < len; i++)
         {
           octave_quit ();
-          result(i) = std::pow (a(i), b(i));
+          result(i) = zzpow (a(i), b(i));
         }
 
       retval = result;
@@ -1304,7 +1255,7 @@ elem_xpow (const NDArray& a, const Compl
   for (octave_idx_type i = 0; i < a.numel (); i++)
     {
       octave_quit ();
-      result(i) = std::pow (a(i), b);
+      result(i) = zzpow (a(i), b);
     }
 
   return result;
@@ -1330,7 +1281,7 @@ elem_xpow (const NDArray& a, const Compl
   for (octave_idx_type i = 0; i < a.numel (); i++)
     {
       octave_quit ();
-      result(i) = std::pow (a(i), b(i));
+      result(i) = zzpow (a(i), b(i));
     }
 
   return result;
@@ -1345,11 +1296,7 @@ elem_xpow (const Complex& a, const NDArr
   for (octave_idx_type i = 0; i < b.numel (); i++)
     {
       octave_quit ();
-      double btmp = b(i);
-      if (xisint (btmp))
-        result(i) = std::pow (a, static_cast<int> (btmp));
-      else
-        result(i) = std::pow (a, btmp);
+      result(i) = zzpow (a, b(i));
     }
 
   return result;
@@ -1364,7 +1311,7 @@ elem_xpow (const Complex& a, const Compl
   for (octave_idx_type i = 0; i < b.numel (); i++)
     {
       octave_quit ();
-      result(i) = std::pow (a, b(i));
+      result(i) = zzpow (a, b(i));
     }
 
   return result;
@@ -1376,29 +1323,17 @@ elem_xpow (const ComplexNDArray& a, doub
 {
   ComplexNDArray result (a.dims ());
 
-  if (xisint (b))
+  if (xisint (b) && static_cast<int> (b) == -1)
     {
-      int bint = static_cast<int> (b);
-      if (bint == -1)
-        {
-          for (octave_idx_type i = 0; i < a.numel (); i++)
-            result.xelem (i) = 1.0 / a(i);
-        }
-      else
-        {
-          for (octave_idx_type i = 0; i < a.numel (); i++)
-            {
-              octave_quit ();
-              result(i) = std::pow (a(i), bint);
-            }
-        }
+      for (octave_idx_type i = 0; i < a.numel (); i++)
+        result.xelem (i) = 1.0 / a(i);
     }
   else
     {
       for (octave_idx_type i = 0; i < a.numel (); i++)
         {
           octave_quit ();
-          result(i) = std::pow (a(i), b);
+          result(i) = zzpow (a(i), b);
         }
     }
 
@@ -1425,11 +1360,7 @@ elem_xpow (const ComplexNDArray& a, cons
   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);
+      result(i) = zzpow (a(i), b(i));
     }
 
   return result;
@@ -1444,7 +1375,7 @@ elem_xpow (const ComplexNDArray& a, cons
   for (octave_idx_type i = 0; i < a.numel (); i++)
     {
       octave_quit ();
-      result(i) = std::pow (a(i), b);
+      result(i) = zzpow (a(i), b);
     }
 
   return result;
@@ -1470,7 +1401,7 @@ elem_xpow (const ComplexNDArray& a, cons
   for (octave_idx_type i = 0; i < a.numel (); i++)
     {
       octave_quit ();
-      result(i) = std::pow (a(i), b(i));
+      result(i) = zzpow (a(i), b(i));
     }
 
   return result;
@@ -1499,10 +1430,10 @@ xpow (float a, float b)
     {
       FloatComplex atmp (a);
 
-      return std::pow (atmp, b);
+      return zzpow (atmp, b);
     }
   else
-    retval = std::pow (a, b);
+    retval = zzpow (a, b);
 
   return retval;
 }
@@ -1530,9 +1461,9 @@ xpow (float a, const FloatMatrix& b)
         {
           FloatComplex elt = lambda(i);
           if (std::imag (elt) == 0.0)
-            lambda(i) = std::pow (a, std::real (elt));
+            lambda(i) = zzpow (a, std::real (elt));
           else
-            lambda(i) = std::pow (a, elt);
+            lambda(i) = zzpow (a, elt);
         }
       FloatComplexDiagMatrix D (lambda);
 
@@ -1555,7 +1486,7 @@ xpow (float a, const FloatMatrix& b)
 octave_value
 xpow (float a, const FloatComplex& b)
 {
-  FloatComplex result = std::pow (a, b);
+  FloatComplex result = zzpow (a, b);
   return result;
 }
 
@@ -1582,9 +1513,9 @@ xpow (float a, const FloatComplexMatrix&
         {
           FloatComplex elt = lambda(i);
           if (std::imag (elt) == 0.0)
-            lambda(i) = std::pow (a, std::real (elt));
+            lambda(i) = zzpow (a, std::real (elt));
           else
-            lambda(i) = std::pow (a, elt);
+            lambda(i) = zzpow (a, elt);
         }
       FloatComplexDiagMatrix D (lambda);
 
@@ -1669,7 +1600,7 @@ xpow (const FloatMatrix& a, float b)
           FloatComplexMatrix Q (a_eig.right_eigenvectors ());
 
           for (octave_idx_type i = 0; i < nr; i++)
-            lambda(i) = std::pow (lambda(i), b);
+            lambda(i) = zzpow (lambda(i), b);
 
           FloatComplexDiagMatrix D (lambda);
 
@@ -1700,14 +1631,14 @@ xpow (const FloatDiagMatrix& a, float b)
     {
       FloatDiagMatrix r (nr, nc);
       for (octave_idx_type i = 0; i < nc; i++)
-        r.dgelem (i) = std::pow (a.dgelem (i), b);
+        r.dgelem (i) = zzpow (a.dgelem (i), b);
       retval = r;
     }
   else
     {
       FloatComplexDiagMatrix r (nr, nc);
       for (octave_idx_type i = 0; i < nc; i++)
-        r.dgelem (i) = std::pow (static_cast<FloatComplex> (a.dgelem (i)), b);
+        r.dgelem (i) = zzpow (static_cast<FloatComplex> (a.dgelem (i)), b);
       retval = r;
     }
 
@@ -1734,7 +1665,7 @@ xpow (const FloatMatrix& a, const FloatC
       FloatComplexMatrix Q (a_eig.right_eigenvectors ());
 
       for (octave_idx_type i = 0; i < nr; i++)
-        lambda(i) = std::pow (lambda(i), b);
+        lambda(i) = zzpow (lambda(i), b);
 
       FloatComplexDiagMatrix D (lambda);
 
@@ -1752,13 +1683,7 @@ xpow (const FloatMatrix& a, const FloatC
 octave_value
 xpow (const FloatComplex& a, float b)
 {
-  FloatComplex result;
-
-  if (xisint (b))
-    result = std::pow (a, static_cast<int> (b));
-  else
-    result = std::pow (a, b);
-
+  FloatComplex result = zzpow (a, b);
   return result;
 }
 
@@ -1785,9 +1710,9 @@ xpow (const FloatComplex& a, const Float
         {
           FloatComplex elt = lambda(i);
           if (std::imag (elt) == 0.0)
-            lambda(i) = std::pow (a, std::real (elt));
+            lambda(i) = zzpow (a, std::real (elt));
           else
-            lambda(i) = std::pow (a, elt);
+            lambda(i) = zzpow (a, elt);
         }
       FloatComplexDiagMatrix D (lambda);
 
@@ -1805,8 +1730,7 @@ xpow (const FloatComplex& a, const Float
 octave_value
 xpow (const FloatComplex& a, const FloatComplex& b)
 {
-  FloatComplex result;
-  result = std::pow (a, b);
+  FloatComplex result = zzpow (a, b);
   return result;
 }
 
@@ -1833,9 +1757,9 @@ xpow (const FloatComplex& a, const Float
         {
           FloatComplex elt = lambda(i);
           if (std::imag (elt) == 0.0)
-            lambda(i) = std::pow (a, std::real (elt));
+            lambda(i) = zzpow (a, std::real (elt));
           else
-            lambda(i) = std::pow (a, elt);
+            lambda(i) = zzpow (a, elt);
         }
       FloatComplexDiagMatrix D (lambda);
 
@@ -1920,7 +1844,7 @@ xpow (const FloatComplexMatrix& a, float
           FloatComplexMatrix Q (a_eig.right_eigenvectors ());
 
           for (octave_idx_type i = 0; i < nr; i++)
-            lambda(i) = std::pow (lambda(i), b);
+            lambda(i) = zzpow (lambda(i), b);
 
           FloatComplexDiagMatrix D (lambda);
 
@@ -1955,7 +1879,7 @@ xpow (const FloatComplexMatrix& a, const
       FloatComplexMatrix Q (a_eig.right_eigenvectors ());
 
       for (octave_idx_type i = 0; i < nr; i++)
-        lambda(i) = std::pow (lambda(i), b);
+        lambda(i) = zzpow (lambda(i), b);
 
       FloatComplexDiagMatrix D (lambda);
 
@@ -1983,7 +1907,7 @@ xpow (const FloatComplexDiagMatrix& a, c
 
   FloatComplexDiagMatrix r (nr, nc);
   for (octave_idx_type i = 0; i < nc; i++)
-    r(i, i) = std::pow (a(i, i), b);
+    r(i, i) = zzpow (a(i, i), b);
   retval = r;
 
   return retval;
@@ -2050,7 +1974,7 @@ elem_xpow (float a, const FloatMatrix& b
         for (octave_idx_type i = 0; i < nr; i++)
           {
             octave_quit ();
-            result(i, j) = std::pow (atmp, b(i, j));
+            result(i, j) = zzpow (atmp, b(i, j));
           }
 
       retval = result;
@@ -2063,7 +1987,7 @@ elem_xpow (float a, const FloatMatrix& b
         for (octave_idx_type i = 0; i < nr; i++)
           {
             octave_quit ();
-            result(i, j) = std::pow (a, b(i, j));
+            result(i, j) = zzpow (a, b(i, j));
           }
 
       retval = result;
@@ -2086,7 +2010,7 @@ elem_xpow (float a, const FloatComplexMa
     for (octave_idx_type i = 0; i < nr; i++)
       {
         octave_quit ();
-        result(i, j) = std::pow (atmp, b(i, j));
+        result(i, j) = zzpow (atmp, b(i, j));
       }
 
   return result;
@@ -2112,7 +2036,7 @@ elem_xpow (const FloatMatrix& a, float b
 
             FloatComplex atmp (a(i, j));
 
-            result(i, j) = std::pow (atmp, b);
+            result(i, j) = zzpow (atmp, b);
           }
 
       retval = result;
@@ -2125,7 +2049,7 @@ elem_xpow (const FloatMatrix& a, float b
         for (octave_idx_type i = 0; i < nr; i++)
           {
             octave_quit ();
-            result(i, j) = std::pow (a(i, j), b);
+            result(i, j) = zzpow (a(i, j), b);
           }
 
       retval = result;
@@ -2175,7 +2099,7 @@ done:
             octave_quit ();
             FloatComplex atmp (a(i, j));
             FloatComplex btmp (b(i, j));
-            complex_result(i, j) = std::pow (atmp, btmp);
+            complex_result(i, j) = zzpow (atmp, btmp);
           }
 
       retval = complex_result;
@@ -2188,7 +2112,7 @@ done:
         for (octave_idx_type i = 0; i < nr; i++)
           {
             octave_quit ();
-            result(i, j) = std::pow (a(i, j), b(i, j));
+            result(i, j) = zzpow (a(i, j), b(i, j));
           }
 
       retval = result;
@@ -2210,7 +2134,7 @@ elem_xpow (const FloatMatrix& a, const F
     for (octave_idx_type i = 0; i < nr; i++)
       {
         octave_quit ();
-        result(i, j) = std::pow (FloatComplex (a(i, j)), b);
+        result(i, j) = zzpow (FloatComplex (a(i, j)), b);
       }
 
   return result;
@@ -2235,7 +2159,7 @@ elem_xpow (const FloatMatrix& a, const F
     for (octave_idx_type i = 0; i < nr; i++)
       {
         octave_quit ();
-        result(i, j) = std::pow (FloatComplex (a(i, j)), b(i, j));
+        result(i, j) = zzpow (FloatComplex (a(i, j)), b(i, j));
       }
 
   return result;
@@ -2254,11 +2178,7 @@ elem_xpow (const FloatComplex& a, const 
     for (octave_idx_type i = 0; i < nr; i++)
       {
         octave_quit ();
-        float btmp = b(i, j);
-        if (xisint (btmp))
-          result(i, j) = std::pow (a, static_cast<int> (btmp));
-        else
-          result(i, j) = std::pow (a, btmp);
+        result(i, j) = zzpow (a, b(i, j));
       }
 
   return result;
@@ -2277,7 +2197,7 @@ elem_xpow (const FloatComplex& a, const 
     for (octave_idx_type i = 0; i < nr; i++)
       {
         octave_quit ();
-        result(i, j) = std::pow (a, b(i, j));
+        result(i, j) = zzpow (a, b(i, j));
       }
 
   return result;
@@ -2292,25 +2212,12 @@ elem_xpow (const FloatComplexMatrix& a, 
 
   FloatComplexMatrix result (nr, nc);
 
-  if (xisint (b))
-    {
-      int bint = static_cast<int> (b);
-      for (octave_idx_type j = 0; j < nc; j++)
-        for (octave_idx_type i = 0; i < nr; i++)
-          {
-            octave_quit ();
-            result(i, j) = std::pow (a(i, j), bint);
-          }
-    }
-  else
-    {
-      for (octave_idx_type j = 0; j < nc; j++)
-        for (octave_idx_type i = 0; i < nr; i++)
-          {
-            octave_quit ();
-            result(i, j) = std::pow (a(i, j), b);
-          }
-    }
+  for (octave_idx_type j = 0; j < nc; j++)
+    for (octave_idx_type i = 0; i < nr; i++)
+      {
+        octave_quit ();
+        result(i, j) = zzpow (a(i, j), b);
+      }
 
   return result;
 }
@@ -2334,11 +2241,7 @@ elem_xpow (const FloatComplexMatrix& a, 
     for (octave_idx_type i = 0; i < nr; i++)
       {
         octave_quit ();
-        float btmp = b(i, j);
-        if (xisint (btmp))
-          result(i, j) = std::pow (a(i, j), static_cast<int> (btmp));
-        else
-          result(i, j) = std::pow (a(i, j), btmp);
+        result(i, j) = zzpow (a(i, j), b(i, j));
       }
 
   return result;
@@ -2357,7 +2260,7 @@ elem_xpow (const FloatComplexMatrix& a, 
     for (octave_idx_type i = 0; i < nr; i++)
       {
         octave_quit ();
-        result(i, j) = std::pow (a(i, j), b);
+        result(i, j) = zzpow (a(i, j), b);
       }
 
   return result;
@@ -2382,7 +2285,7 @@ elem_xpow (const FloatComplexMatrix& a, 
     for (octave_idx_type i = 0; i < nr; i++)
       {
         octave_quit ();
-        result(i, j) = std::pow (a(i, j), b(i, j));
+        result(i, j) = zzpow (a(i, j), b(i, j));
       }
 
   return result;
@@ -2429,7 +2332,7 @@ elem_xpow (float a, const FloatNDArray& 
       for (octave_idx_type i = 0; i < b.numel (); i++)
         {
           octave_quit ();
-          result(i) = std::pow (atmp, b(i));
+          result(i) = zzpow (atmp, b(i));
         }
 
       retval = result;
@@ -2440,7 +2343,7 @@ elem_xpow (float a, const FloatNDArray& 
       for (octave_idx_type i = 0; i < b.numel (); i++)
         {
           octave_quit ();
-          result(i) = std::pow (a, b(i));
+          result(i) = zzpow (a, b(i));
         }
 
       retval = result;
@@ -2458,7 +2361,7 @@ elem_xpow (float a, const FloatComplexND
   for (octave_idx_type i = 0; i < b.numel (); i++)
     {
       octave_quit ();
-      result(i) = std::pow (a, b(i));
+      result(i) = zzpow (a, b(i));
     }
 
   return result;
@@ -2482,7 +2385,7 @@ elem_xpow (const FloatNDArray& a, float 
 
               FloatComplex atmp (a(i));
 
-              result(i) = std::pow (atmp, b);
+              result(i) = zzpow (atmp, b);
             }
 
           retval = result;
@@ -2493,7 +2396,7 @@ elem_xpow (const FloatNDArray& a, float 
           for (octave_idx_type i = 0; i < a.numel (); i++)
             {
               octave_quit ();
-              result(i) = std::pow (a(i), b);
+              result(i) = zzpow (a(i), b);
             }
 
           retval = result;
@@ -2524,7 +2427,7 @@ elem_xpow (const FloatNDArray& a, float 
           for (octave_idx_type i = 0; i < a.numel (); i++)
             {
               octave_quit ();
-              result(i) = std::pow (a(i), ib);
+              result(i) = zzpow (a(i), ib);
             }
         }
 
@@ -2583,7 +2486,7 @@ done:
         {
           octave_quit ();
           FloatComplex atmp (a(i));
-          complex_result(i) = std::pow (atmp, b(i));
+          complex_result(i) = zzpow (atmp, b(i));
         }
 
       retval = complex_result;
@@ -2595,7 +2498,7 @@ done:
       for (octave_idx_type i = 0; i < len; i++)
         {
           octave_quit ();
-          result(i) = std::pow (a(i), b(i));
+          result(i) = zzpow (a(i), b(i));
         }
 
       retval = result;
@@ -2613,7 +2516,7 @@ elem_xpow (const FloatNDArray& a, const 
   for (octave_idx_type i = 0; i < a.numel (); i++)
     {
       octave_quit ();
-      result(i) = std::pow (a(i), b);
+      result(i) = zzpow (a(i), b);
     }
 
   return result;
@@ -2639,7 +2542,7 @@ elem_xpow (const FloatNDArray& a, const 
   for (octave_idx_type i = 0; i < a.numel (); i++)
     {
       octave_quit ();
-      result(i) = std::pow (a(i), b(i));
+      result(i) = zzpow (a(i), b(i));
     }
 
   return result;
@@ -2654,11 +2557,7 @@ elem_xpow (const FloatComplex& a, const 
   for (octave_idx_type i = 0; i < b.numel (); i++)
     {
       octave_quit ();
-      float btmp = b(i);
-      if (xisint (btmp))
-        result(i) = std::pow (a, static_cast<int> (btmp));
-      else
-        result(i) = std::pow (a, btmp);
+      result(i) = zzpow (a, b(i));
     }
 
   return result;
@@ -2673,7 +2572,7 @@ elem_xpow (const FloatComplex& a, const 
   for (octave_idx_type i = 0; i < b.numel (); i++)
     {
       octave_quit ();
-      result(i) = std::pow (a, b(i));
+      result(i) = zzpow (a, b(i));
     }
 
   return result;
@@ -2685,29 +2584,17 @@ elem_xpow (const FloatComplexNDArray& a,
 {
   FloatComplexNDArray result (a.dims ());
 
-  if (xisint (b))
+  if (xisint (b) && static_cast<int> (b) == -1)
     {
-      int bint = static_cast<int> (b);
-      if (bint == -1)
-        {
-          for (octave_idx_type i = 0; i < a.numel (); i++)
-            result.xelem (i) = 1.0f / a(i);
-        }
-      else
-        {
-          for (octave_idx_type i = 0; i < a.numel (); i++)
-            {
-              octave_quit ();
-              result(i) = std::pow (a(i), bint);
-            }
-        }
+      for (octave_idx_type i = 0; i < a.numel (); i++)
+        result.xelem (i) = 1.0f / a(i);
     }
   else
     {
       for (octave_idx_type i = 0; i < a.numel (); i++)
         {
           octave_quit ();
-          result(i) = std::pow (a(i), b);
+          result(i) = zzpow (a(i), b);
         }
     }
 
@@ -2734,11 +2621,7 @@ elem_xpow (const FloatComplexNDArray& a,
   for (octave_idx_type i = 0; i < a.numel (); i++)
     {
       octave_quit ();
-      float btmp = b(i);
-      if (xisint (btmp))
-        result(i) = std::pow (a(i), static_cast<int> (btmp));
-      else
-        result(i) = std::pow (a(i), btmp);
+      result(i) = zzpow (a(i), b(i));
     }
 
   return result;
@@ -2753,7 +2636,7 @@ elem_xpow (const FloatComplexNDArray& a,
   for (octave_idx_type i = 0; i < a.numel (); i++)
     {
       octave_quit ();
-      result(i) = std::pow (a(i), b);
+      result(i) = zzpow (a(i), b);
     }
 
   return result;
@@ -2779,7 +2662,7 @@ elem_xpow (const FloatComplexNDArray& a,
   for (octave_idx_type i = 0; i < a.numel (); i++)
     {
       octave_quit ();
-      result(i) = std::pow (a(i), b(i));
+      result(i) = zzpow (a(i), b(i));
     }
 
   return result;
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,38 +404,91 @@ 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.
+// The C++ std::pow(0,0) and pow(0,0) calls do not necessarily return a
+// value of 1.  Force that to be the case with an intermediate function.
+
+template<class T> std::complex<T>
+inline zzpow (const std::complex<T>& x, const std::complex<T>& y)
+{
+  if (y.real () == 0 && y.imag () == 0)
+    return std::complex<T> (1, 0);
+  else
+    return std::pow (x, y);
+}
+
+template<class T> std::complex<T>
+inline zzpow (const std::complex<T>& x, const T& y)
+{
+  if (y == 0)
+    return std::complex<T> (1, 0);
+  else
+    return std::pow (x, y);
+}
+
+template<class T> std::complex<T>
+inline zzpow (const T& x, const std::complex<T>& y)
+{
+  if (y.real () == 0 && y.imag () == 0)
+    return std::complex<T> (1, 0);
+  else
+    return std::pow (x, y);
+}
+
+double inline zzpow (double base, double exponent)
+{
+  if (exponent == 0.0)
+    return 1.0;
+  else
+    return pow (base, exponent);
+}
+
+float inline zzpow (float base, float exponent)
+{
+  if (exponent == 0.0)
+    return 1.0;
+  else
+    return pow (base, exponent);
+}
+
+long double inline zzpow (long double base, long double exponent)
+{
+  if (exponent == 0.0L)
+    return 1.0;
+  else
+    return pow (base, exponent);
+}
+
+template <typename Type1, typename Type2>
+double inline zzpow (Type1 base, Type2 exponent)
+{
+  if (exponent == Type2 (0))
+    return 1.0;
+  else
+    return pow (base, exponent);
+}
 
 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]);
+    r[i] = zzpow (x[i], y[i]);
 }
 
 template <typename R, typename X, typename Y>
 inline void
 mx_inline_pow (size_t n, R *r, const X *x, Y y)
 {
-  using std::pow;
-
   for (size_t i = 0; i < n; i++)
-    r[i] = pow (x[i], y);
+    r[i] = zzpow (x[i], y);
 }
 
 template <typename R, typename X, typename Y>
 inline void
 mx_inline_pow (size_t n, R *r, X x, const Y *y)
 {
-  using std::pow;
-
   for (size_t i = 0; i < n; i++)
-    r[i] = pow (x, y[i]);
+    r[i] = zzpow (x, y[i]);
 }
 
 // Arbitrary function appliers.
