Mercurial > hg > octave-nkf
diff src/xpow.cc @ 9103:10bed8fbec99
optimize scalar .^ range operation
author | Jaroslav Hajek <highegg@gmail.com> |
---|---|
date | Wed, 08 Apr 2009 10:53:06 +0200 |
parents | a7c00773a089 |
children | f39b98237d5c |
line wrap: on
line diff
--- 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)