Mercurial > hg > octave-lojdl
changeset 9123:f39b98237d5c
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 +++ b/src/ChangeLog @@ -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 +++ b/src/xpow.cc @@ -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; }