Mercurial > octave

changeset 29816:f54cfd60725e

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Fix mpower for additional cases where negative scalar base is used (bug #60786). * xpow.cc (xpow (double, ComplexMatrix), xpow (float, FloatComplexMatrix)): 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 Thu, 24 Jun 2021 13:51:52 -0700
parents e637c2342433
children 0511e3638724
files libinterp/corefcn/xpow.cc
diffstat 1 files changed, 42 insertions(+), 12 deletions(-) [+]
line wrap: on
line diff
--- a/libinterp/corefcn/xpow.cc	Thu Jun 24 20:47:08 2021 +0200
+++ b/libinterp/corefcn/xpow.cc	Thu Jun 24 13:51:52 2021 -0700
@@ -201,13 +201,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);
 
@@ -1650,13 +1665,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);