changeset 9103:10bed8fbec99

optimize scalar .^ range operation
author Jaroslav Hajek <highegg@gmail.com>
date Wed, 08 Apr 2009 10:53:06 +0200
parents 9dc516d36175
children e0250e2b60ed
files src/ChangeLog src/OPERATORS/op-range.cc src/xpow.cc src/xpow.h
diffstat 4 files changed, 79 insertions(+), 0 deletions(-) [+]
line wrap: on
line diff
--- a/src/ChangeLog
+++ b/src/ChangeLog
@@ -1,3 +1,11 @@
+2009-04-08  Jaroslav Hajek  <highegg@gmail.com>
+
+	* xpow.cc (elem_xpow (double, const Range&),
+	elem_xpow (const Complex&, const Range&)): New functions.
+	* xpow.h: Declare them.
+	* OPERATORS/op-range.cc: Define scalar .^ range and complex .^ range &
+	install them.
+
 2009-04-04  Jaroslav Hajek  <highegg@gmail.com>
 
 	* ov-struct.cc (octave_struct::subsasgn): Fix reference counting
--- a/src/OPERATORS/op-range.cc
+++ b/src/OPERATORS/op-range.cc
@@ -40,6 +40,7 @@
 #include "ov-typeinfo.h"
 #include "ov-null-mat.h"
 #include "ops.h"
+#include "xpow.h"
 
 // range unary ops.
 
@@ -67,6 +68,9 @@
 DEFBINOP_OP (mulrs, range, scalar, *)
 DEFBINOP_OP (mulsr, scalar, range, *)
 
+DEFBINOP_FN (el_powsr, scalar, range, elem_xpow)
+DEFBINOP_FN (el_powcsr, complex, range, elem_xpow)
+
 DEFNDCATOP_FN (r_r, range, range, array, array, concat)
 DEFNDCATOP_FN (r_s, range, scalar, array, array, concat)
 DEFNDCATOP_FN (r_m, range, matrix, array, array, concat)
@@ -106,6 +110,9 @@
   INSTALL_BINOP (op_mul, octave_range, octave_scalar, mulrs);
   INSTALL_BINOP (op_mul, octave_scalar, octave_range, mulsr);
 
+  INSTALL_BINOP (op_el_pow, octave_scalar, octave_range, el_powsr);
+  INSTALL_BINOP (op_el_pow, octave_complex, octave_range, el_powcsr);
+
   INSTALL_CATOP (octave_range, octave_range, r_r);
   INSTALL_CATOP (octave_range, octave_scalar, r_s);
   INSTALL_CATOP (octave_range, octave_matrix, r_m);
--- a/src/xpow.cc
+++ b/src/xpow.cc
@@ -41,6 +41,7 @@
 #include "PermMatrix.h"
 #include "mx-cm-cdm.h"
 #include "oct-cmplx.h"
+#include "Range.h"
 #include "quit.h"
 
 #include "error.h"
@@ -691,6 +692,43 @@
   return result;
 }
 
+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 ())
+    {
+      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;
+    }
+  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;
+    }
+
+  return retval;
+}
+
 // -*- 3 -*-
 octave_value
 elem_xpow (const Matrix& a, double b)
@@ -888,6 +926,29 @@
   return result;
 }
 
+octave_value
+elem_xpow (const Complex& a, const Range& r)
+{
+  octave_value retval;
+
+  if (r.nelem () <= 0)
+    retval = Matrix ();
+  else
+    {
+      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;
+    }
+
+  return retval;
+}
+
 // -*- 9 -*-
 octave_value
 elem_xpow (const ComplexMatrix& a, double b)
--- a/src/xpow.h
+++ b/src/xpow.h
@@ -40,6 +40,7 @@
 class ComplexNDArray;
 class FloatComplexNDArray;
 class octave_value;
+class Range;
 
 extern octave_value xpow (double a, double b);
 extern octave_value xpow (double a, const Matrix& b);
@@ -67,6 +68,7 @@
 
 extern octave_value elem_xpow (double a, const Matrix& b);
 extern octave_value elem_xpow (double a, const ComplexMatrix& b);
+extern octave_value elem_xpow (double a, const Range& r);
 
 extern octave_value elem_xpow (const Matrix& a, double b);
 extern octave_value elem_xpow (const Matrix& a, const Matrix& b);
@@ -75,6 +77,7 @@
 
 extern octave_value elem_xpow (const Complex& a, const Matrix& b);
 extern octave_value elem_xpow (const Complex& a, const ComplexMatrix& b);
+extern octave_value elem_xpow (const Complex& a, const Range& r);
 
 extern octave_value elem_xpow (const ComplexMatrix& a, double b);
 extern octave_value elem_xpow (const ComplexMatrix& a, const Matrix& b);