changeset 29803:540f25090412

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.
author Rik <rik@octave.org>
date Tue, 22 Jun 2021 13:53:49 -0700
parents 10408a66305b
children 977e4ad8966d
files libinterp/corefcn/xpow.cc
diffstat 1 files changed, 42 insertions(+), 12 deletions(-) [+]
line wrap: on
line diff
--- 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);