The same patch for GiNaC 1.3. Make code converting terms into primitive polynomials (i.e. (-x+y)*(x-y) -> -(x-y)^2, (3*x+y)/(6*x+2*y)^2 -> 1/4*(y+3*x)^(-1)) a little bit more efficient: the numeric coefficient is calculated in one sweep, sums divieded by numbers in place. --- ginac/compiler.h | 12 ++++++++++++ ginac/mul.cpp | 48 ++++++++++++++++++++++++++++++++++-------------- ginac/power.cpp | 34 ++++++++++++++++++++-------------- 3 files changed, 66 insertions(+), 28 deletions(-) diff --git a/ginac/compiler.h b/ginac/compiler.h new file mode 100644 index 0000000..66f6e92 --- /dev/null +++ b/ginac/compiler.h @@ -0,0 +1,12 @@ +#ifndef GINAC_COMPILER_DEP_HH +#define GINAC_COMPILER_DEP_HH + +#ifdef __GNUC__ +#define unlikely(cond) __builtin_expect((cond), 0) +#define likely(cond) __builtin_expect((cond), 1) +#else +#define unlikely(cond) (cond) +#define likely(cond) (cond) +#endif + +#endif /* GINAC_COMPILER_DEP_HH */ diff --git a/ginac/mul.cpp b/ginac/mul.cpp index 230717a..3c87d5a 100644 --- a/ginac/mul.cpp +++ b/ginac/mul.cpp @@ -34,6 +34,7 @@ #include "lst.h" #include "archive.h" #include "utils.h" +#include "compiler.h" namespace GiNaC { @@ -414,7 +415,7 @@ ex mul::eval(int level) const return *this; } - int seq_size = seq.size(); + size_t seq_size = seq.size(); if (overall_coeff.is_zero()) { // *(...,x;0) -> 0 return _ex0; @@ -440,14 +441,18 @@ ex mul::eval(int level) const ex_to<numeric>(addref.overall_coeff). mul_dyn(ex_to<numeric>(overall_coeff)))) ->setflag(status_flags::dynallocated | status_flags::evaluated); - } else if (seq_size >= 2) { + } else if ((seq_size >= 2) && (! (flags & status_flags::expanded))) { // Strip the content and the unit part from each term. Thus // things like (-x+a)*(3*x-3*a) automagically turn into - 3*(x-a)2 epvector::const_iterator last = seq.end(); epvector::const_iterator i = seq.begin(); + epvector::const_iterator j = seq.begin(); + std::auto_ptr<epvector> s(new epvector); + numeric oc = *_num1_p; + bool something_changed = false; while (i!=last) { - if (! (is_a<add>(i->rest) && i->coeff.is_equal(_ex1))) { + if (likely(! (is_a<add>(i->rest) && i->coeff.is_equal(_ex1)))) { // power::eval has such a rule, no need to handle powers here ++i; continue; @@ -464,32 +469,47 @@ ex mul::eval(int level) const // very unlucky event it can even loop forever). Hopefully the main // variable will be the same for all terms in *this const bool unit_normal = lead_coeff.is_pos_integer(); - if ((c == *_num1_p) && ((! canonicalizable) || unit_normal)) { + if (likely((c == *_num1_p) && ((! canonicalizable) || unit_normal))) { ++i; continue; } - std::auto_ptr<epvector> s(new epvector); - s->reserve(seq.size()); + if (! something_changed) { + s->reserve(seq_size); + something_changed = true; + } - epvector::const_iterator j=seq.begin(); - while (j!=i) { + while ((j!=i) && (j!=last)) { s->push_back(*j); ++j; } - if (! unit_normal) { + if (! unit_normal) c = c.mul(*_num_1_p); - } - const ex primitive = (i->rest)/c; - s->push_back(expair(primitive, _ex1)); - ++j; + oc = oc.mul(c); + + // divide add by the number in place to save at least 2 .eval() calls + const add& addref = ex_to<add>(i->rest); + add* primitive = new add(addref); + primitive->setflag(status_flags::dynallocated); + primitive->clearflag(status_flags::hash_calculated); + primitive->overall_coeff = ex_to<numeric>(primitive->overall_coeff).div_dyn(c); + for (epvector::iterator ai = primitive->seq.begin(); + ai != primitive->seq.end(); ++ai) + ai->coeff = ex_to<numeric>(ai->coeff).div_dyn(c); + + s->push_back(expair(*primitive, _ex1)); + + ++i; + ++j; + } + if (something_changed) { while (j!=last) { s->push_back(*j); ++j; } - return (new mul(s, ex_to<numeric>(overall_coeff).mul_dyn(c)) + return (new mul(s, ex_to<numeric>(overall_coeff).mul_dyn(oc)) )->setflag(status_flags::dynallocated); } } diff --git a/ginac/power.cpp b/ginac/power.cpp index 0dfbc8f..a94cc47 100644 --- a/ginac/power.cpp +++ b/ginac/power.cpp @@ -40,6 +40,8 @@ #include "lst.h" #include "archive.h" #include "utils.h" +#include "relational.h" +#include "compiler.h" namespace GiNaC { @@ -465,25 +467,29 @@ ex power::eval(int level) const // (2*x + 6*y)^(-4) -> 1/16*(x + 3*y)^(-4) if (num_exponent->is_integer() && is_exactly_a<add>(ebasis)) { - const numeric icont = ebasis.integer_content(); + numeric icont = ebasis.integer_content(); const numeric& lead_coeff = ex_to<numeric>(ex_to<add>(ebasis).seq.begin()->coeff).div_dyn(icont); const bool canonicalizable = lead_coeff.is_integer(); const bool unit_normal = lead_coeff.is_pos_integer(); - - if (icont != *_num1_p) { - return (new mul(power(ebasis/icont, *num_exponent), power(icont, *num_exponent)) - )->setflag(status_flags::dynallocated); - } - - if (canonicalizable && (! unit_normal)) { - if (num_exponent->is_even()) { - return power(-ebasis, *num_exponent); - } else { - return (new mul(power(-ebasis, *num_exponent), *_num_1_p) - )->setflag(status_flags::dynallocated); - } + if (canonicalizable && (! unit_normal)) + icont = icont.mul(*_num_1_p); + + if (canonicalizable && (icont != *_num1_p)) { + const add& addref = ex_to<add>(ebasis); + add* addp = new add(addref); + addp->setflag(status_flags::dynallocated); + addp->clearflag(status_flags::hash_calculated); + addp->overall_coeff = ex_to<numeric>(addp->overall_coeff).div_dyn(icont); + for (epvector::iterator i = addp->seq.begin(); i != addp->seq.end(); ++i) + i->coeff = ex_to<numeric>(i->coeff).div_dyn(icont); + + const numeric c = icont.power(*num_exponent); + if (likely(c != *_num1_p)) + return (new mul(power(*addp, *num_exponent), c))->setflag(status_flags::dynallocated); + else + return power(*addp, *num_exponent); } } -- 1.4.4.4 Best regards, Alexei -- All science is either physics or stamp collecting.