Mercurial > octave-nkf

changeset 9123:f39b98237d5c

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
use more robust and less aggressive scalar .^ range optimization
author Jaroslav Hajek <highegg@gmail.com>
date Thu, 16 Apr 2009 07:02:31 +0200
parents 8ca06fd9c6ef
children 47f19c11b558
files src/ChangeLog src/xpow.cc
diffstat 2 files changed, 70 insertions(+), 39 deletions(-) [+]
line wrap: on
line diff
--- a/src/ChangeLog	Wed Apr 15 12:29:31 2009 +0200
+++ b/src/ChangeLog	Thu Apr 16 07:02:31 2009 +0200
@@ -1,3 +1,10 @@
+2009-04-16  Jaroslav Hajek  <highegg@gmail.com>
+
+	* xpow.cc (same_sign): New helper function.
+	(elem_xpow (double, const Range&), elem_xpow (const Complex&, const Range&)): 
+	Only optimize monotonic-magnitude integer ranges, start always from
+	the smaller end.
+
 2008-04-11  David Bateman  <dbateman@free.fr>
 
 	* ov-cell.cc (Fstruct2cell): Treat possible trailing singleton for
--- a/src/xpow.cc	Wed Apr 15 12:29:31 2009 +0200
+++ b/src/xpow.cc	Thu Apr 16 07:02:31 2009 +0200
@@ -692,39 +692,46 @@
   return result;
 }
 
+static inline bool 
+same_sign (double a, double b)
+{
+  return (a >= 0 && b >= 0) || (a <= 0 && b <= 0);
+}
+
 octave_value
 elem_xpow (double a, const Range& r)
 {
   octave_value retval;
 
-  if (r.nelem () <= 0)
-    retval = Matrix ();
-  else if (a < 0.0 && ! r.all_elements_are_ints ())
+  // Only optimize powers with ranges that are integer and monotonic in
+  // magnitude.
+  if (r.nelem () > 1 && r.all_elements_are_ints ()
+      && same_sign (r.base (), r.limit ()))
     {
-      ComplexMatrix mat (1, r.nelem ());
-      Complex atmp (a);
-      Complex *pmat = mat.fortran_vec ();
-
-      pmat[0] = std::pow (atmp, r.base ());
-      Complex mul = std::pow (atmp, r.inc ());
-      for (octave_idx_type i = 1; i < r.nelem (); i++)
-        pmat[i] = pmat[i-1] * mul;
-
-      retval = mat;
-    }
+      octave_idx_type n = r.nelem ();
+      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 ());
+          result(0) = base;
+          for (octave_idx_type i = 1; i < n; i++)
+            result(i) = (base *= inc);
+        }
+      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 ());
+          result(n-1) = limit;
+          for (octave_idx_type i = n-2; i >= 0; i--)
+            result(i) = (limit *= inc);
+        }
+
+      retval = result;
+    }  
   else
-    {
-      Matrix mat (1, r.nelem ());
-      double *pmat = mat.fortran_vec ();
-
-      double base = std::pow (a, r.base ());
-      pmat[0] = base;
-      double mul = std::pow (a, r.inc ());
-      for (octave_idx_type i = 1; i < r.nelem (); i++)
-        pmat[i] = (base *= mul);
-
-      retval = mat;
-    }
+    retval = elem_xpow (a, r.matrix_value ());
 
   return retval;
 }
@@ -931,20 +938,37 @@
 {
   octave_value retval;
 
-  if (r.nelem () <= 0)
-    retval = Matrix ();
-  else
+  // Only optimize powers with ranges that are integer and monotonic in
+  // magnitude.
+  if (r.nelem () > 1 && r.all_elements_are_ints ()
+      && same_sign (r.base (), r.limit ()))
     {
-      ComplexMatrix mat (1, r.nelem ());
-      Complex *pmat = mat.fortran_vec ();
-
-      pmat[0] = std::pow (a, r.base ());
-      Complex mul = std::pow (a, r.inc ());
-      for (octave_idx_type i = 1; i < r.nelem (); i++)
-        pmat[i] = pmat[i-1] * mul;
-
-      retval = mat;
-    }
+      octave_idx_type n = r.nelem ();
+      ComplexMatrix result (1, n);
+
+      if (same_sign (r.base (), r.inc ()))
+        {
+          Complex base = std::pow (a, r.base ());
+          Complex inc = std::pow (a, r.inc ());
+          result(0) = base;
+          for (octave_idx_type i = 1; i < n; i++)
+            result(i) = (base *= inc);
+        }
+      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 ());
+          result(n-1) = limit;
+          for (octave_idx_type i = n-2; i >= 0; i--)
+            result(i) = (limit *= inc);
+        }
+
+      retval = result;
+    }  
+  else
+    retval = elem_xpow (a, r.matrix_value ());
+
 
   return retval;
 }