# HG changeset patch
# User Rik <rik@octave.org>
# Date 1624395229 25200
# Node ID 540f25090412f759c1b17516319e8c14eb6bea33
# Parent  10408a66305b71d8cc1311a5c650a7c7eab7a9d8
Fix mpower when negative scalar base is used (bug #60786)

* xpow.cc (xpow (double, Matrix), xpow (float, FloatMatrix)):
Check for base < 0.0 and if found use std::pow with std::complex first argument
so that math library chooses a different template for function which can handle
negative bases.

diff -r 10408a66305b -r 540f25090412 libinterp/corefcn/xpow.cc
--- a/libinterp/corefcn/xpow.cc	Tue Jun 22 11:08:31 2021 -0400
+++ b/libinterp/corefcn/xpow.cc	Tue Jun 22 13:53:49 2021 -0700
@@ -132,13 +132,28 @@
       ComplexColumnVector lambda (b_eig.eigenvalues ());
       ComplexMatrix Q (b_eig.right_eigenvectors ());
 
-      for (octave_idx_type i = 0; i < nr; i++)
+      if (a < 0.0)
         {
-          Complex elt = lambda(i);
-          if (std::imag (elt) == 0.0)
-            lambda(i) = std::pow (a, std::real (elt));
-          else
-            lambda(i) = std::pow (a, elt);
+          Complex acplx (a);
+          for (octave_idx_type i = 0; i < nr; i++)
+            {
+              Complex elt = lambda(i);
+              if (std::imag (elt) == 0.0)
+                lambda(i) = std::pow (acplx, std::real (elt));
+              else
+                lambda(i) = std::pow (a, elt);
+            }
+        }
+      else
+        {
+          for (octave_idx_type i = 0; i < nr; i++)
+            {
+              Complex elt = lambda(i);
+              if (std::imag (elt) == 0.0)
+                lambda(i) = std::pow (a, std::real (elt));
+              else
+                lambda(i) = std::pow (a, elt);
+            }
         }
       ComplexDiagMatrix D (lambda);
 
@@ -1565,13 +1580,28 @@
       FloatComplexColumnVector lambda (b_eig.eigenvalues ());
       FloatComplexMatrix Q (b_eig.right_eigenvectors ());
 
-      for (octave_idx_type i = 0; i < nr; i++)
+      if (a < 0.0)
         {
-          FloatComplex elt = lambda(i);
-          if (std::imag (elt) == 0.0)
-            lambda(i) = std::pow (a, std::real (elt));
-          else
-            lambda(i) = std::pow (a, elt);
+          FloatComplex acplx (a);
+          for (octave_idx_type i = 0; i < nr; i++)
+            {
+              FloatComplex elt = lambda(i);
+              if (std::imag (elt) == 0.0)
+                lambda(i) = std::pow (acplx, std::real (elt));
+              else
+                lambda(i) = std::pow (a, elt);
+            }
+        }
+      else
+        {
+          for (octave_idx_type i = 0; i < nr; i++)
+            {
+              FloatComplex elt = lambda(i);
+              if (std::imag (elt) == 0.0)
+                lambda(i) = std::pow (a, std::real (elt));
+              else
+                lambda(i) = std::pow (a, elt);
+            }
         }
       FloatComplexDiagMatrix D (lambda);