Restore the speed of polynomial operations.
Hello! These patches commit 0633ed8082961673eedc092689e06fa39d6bc322 Author: Jens Vollinga <vollinga@thep.physik.uni-mainz.de> Date: Thu May 19 18:10:40 2005 +0000 Fixed bug in expanding expressions containing dummy indices. [V.Kisil] (see <http://www.ginac.de/pipermail/ginac-cvs/2005-May/001586.html>) commit 6b1b11264fe4872132fef8cb8b22cb983ba01736 Author: Chris Dams <Chris.Dams@mi.infn.it> Date: Thu Jul 13 20:02:38 2006 +0000 Needed to call get_all_dummy_indices_safely more often. (see http://www.ginac.de/pipermail/ginac-cvs/2006-July/001780.html) commit 4cf53676d7823ec3c55a90df7723f7f1492a95b2 Author: Chris Dams <Chris.Dams@mi.infn.it> Date: Mon Jul 24 16:47:43 2006 +0000 Improved dummy index renaming. (see http://www.ginac.de/pipermail/ginac-cvs/2006-July/001782.html) made polynomial operations (in particular expand()) substantially slower. For example, Fateman's polynomial expanding benchark degraded from 32sec. to 80sec. My own code which uses multivariate GCD (which in turn uses .expand() extensively) become 4x slower. I'll post 4 patches which address the issue in reply to this message. Comments, questions, suggestions, etc. are wellcome. Best regards, Alexei -- All science is either physics or stamp collecting.
info(info_flags::has_indices) now works for sums and products. It returns true if the expression has indices (no matter dummy or free), and false otherwise. --- ginac/expairseq.cpp | 23 +++++++++++++++++++++-- ginac/flags.h | 4 +++- ginac/idx.cpp | 7 +++++-- ginac/power.cpp | 15 +++++++++++++++ ginac/symbol.cpp | 35 +++++++++++++++++++---------------- 5 files changed, 63 insertions(+), 21 deletions(-) diff --git a/ginac/expairseq.cpp b/ginac/expairseq.cpp index 757f3c3..79c741a 100644 --- a/ginac/expairseq.cpp +++ b/ginac/expairseq.cpp @@ -24,6 +24,7 @@ #include <algorithm> #include <string> #include <stdexcept> +#include <iterator> #include "expairseq.h" #include "lst.h" @@ -267,8 +268,26 @@ void expairseq::do_print_tree(const print_tree & c, unsigned level) const bool expairseq::info(unsigned inf) const { - if (inf == info_flags::expanded) - return (flags & status_flags::expanded); + switch(inf) { + case info_flags::expanded: + return (flags & status_flags::expanded); + case info_flags::has_indices: { + if (flags & status_flags::has_indices) + return true; + else if (flags & status_flags::has_no_indices) + return false; + for (epvector::const_iterator i = seq.begin(); i != seq.end(); ++i) { + if (i->rest.info(info_flags::has_indices)) { + this->setflag(status_flags::has_indices); + this->clearflag(status_flags::has_no_indices); + return true; + } + } + this->clearflag(status_flags::has_indices); + this->setflag(status_flags::has_no_indices); + return false; + } + } return inherited::info(inf); } diff --git a/ginac/flags.h b/ginac/flags.h index 86d8238..372b899 100644 --- a/ginac/flags.h +++ b/ginac/flags.h @@ -193,7 +193,9 @@ public: evaluated = 0x0002, ///< .eval() has already done its job expanded = 0x0004, ///< .expand(0) has already done its job (other expand() options ignore this flag) hash_calculated = 0x0008, ///< .calchash() has already done its job - not_shareable = 0x0010 ///< don't share instances of this object between different expressions unless explicitly asked to (used by ex::compare()) + not_shareable = 0x0010, ///< don't share instances of this object between different expressions unless explicitly asked to (used by ex::compare()) + has_indices = 0x0020, + has_no_indices = 0x0040 // ! (has_indices || has_no_indices) means "don't know" }; }; diff --git a/ginac/idx.cpp b/ginac/idx.cpp index 3e17569..e966b14 100644 --- a/ginac/idx.cpp +++ b/ginac/idx.cpp @@ -234,8 +234,11 @@ void spinidx::do_print_tree(const print_tree & c, unsigned level) const bool idx::info(unsigned inf) const { - if (inf == info_flags::idx) - return true; + switch(inf) { + case info_flags::idx: + case info_flags::has_indices: + return true; + } return inherited::info(inf); } diff --git a/ginac/power.cpp b/ginac/power.cpp index b490b8c..fe7132a 100644 --- a/ginac/power.cpp +++ b/ginac/power.cpp @@ -240,6 +240,21 @@ bool power::info(unsigned inf) const basis.info(inf); case info_flags::expanded: return (flags & status_flags::expanded); + case info_flags::has_indices: { + if (flags & status_flags::has_indices) + return true; + else if (flags & status_flags::has_no_indices) + return false; + else if (basis.info(info_flags::has_indices)) { + setflag(status_flags::has_indices); + clearflag(status_flags::has_no_indices); + return true; + } else { + clearflag(status_flags::has_indices); + setflag(status_flags::has_no_indices); + return false; + } + } } return inherited::info(inf); } diff --git a/ginac/symbol.cpp b/ginac/symbol.cpp index 9f3b1c6..b68fed8 100644 --- a/ginac/symbol.cpp +++ b/ginac/symbol.cpp @@ -207,22 +207,25 @@ void symbol::do_print_python_repr(const print_python_repr & c, unsigned level) c bool symbol::info(unsigned inf) const { - if (inf == info_flags::symbol) - return true; - if (inf == info_flags::polynomial || - inf == info_flags::integer_polynomial || - inf == info_flags::cinteger_polynomial || - inf == info_flags::rational_polynomial || - inf == info_flags::crational_polynomial || - inf == info_flags::rational_function || - inf == info_flags::expanded) - return true; - if (inf == info_flags::real) - return domain==domain::real || domain==domain::positive; - if (inf == info_flags::positive || inf == info_flags::nonnegative) - return domain == domain::positive; - else - return inherited::info(inf); + switch (inf) { + case info_flags::symbol: + case info_flags::polynomial: + case info_flags::integer_polynomial: + case info_flags::cinteger_polynomial: + case info_flags::rational_polynomial: + case info_flags::crational_polynomial: + case info_flags::rational_function: + case info_flags::expanded: + return true; + case info_flags::real: + return domain == domain::real || domain == domain::positive; + case info_flags::positive: + case info_flags::nonnegative: + return domain == domain::positive; + case info_flags::has_indices: + return false; + } + return inherited::info(inf); } ex symbol::eval(int level) const -- 1.5.2.4 Best regards, Alexei -- All science is either physics or stamp collecting.
Patially solves performance regression in expand(), gcd(), etc. --- ginac/mul.cpp | 56 ++++++++++++++++++++++++++++++++++---------------------- 1 files changed, 34 insertions(+), 22 deletions(-) diff --git a/ginac/mul.cpp b/ginac/mul.cpp index 9ed27c6..db8b9f1 100644 --- a/ginac/mul.cpp +++ b/ginac/mul.cpp @@ -992,6 +992,7 @@ bool mul::can_be_further_expanded(const ex & e) ex mul::expand(unsigned options) const { + const bool skip_idx_rename = ! info(info_flags::has_indices); // First, expand the children std::auto_ptr<epvector> expanded_seqp = expandchildren(options); const epvector & expanded_seq = (expanded_seqp.get() ? *expanded_seqp : seq); @@ -1047,28 +1048,34 @@ ex mul::expand(unsigned options) const ex tmp_accu = (new add(distrseq, add1.overall_coeff*add2.overall_coeff))->setflag(status_flags::dynallocated); exvector add1_dummy_indices, add2_dummy_indices, add_indices; + lst dummy_subs; - for (epvector::const_iterator i=add1begin; i!=add1end; ++i) { - add_indices = get_all_dummy_indices_safely(i->rest); - add1_dummy_indices.insert(add1_dummy_indices.end(), add_indices.begin(), add_indices.end()); - } - for (epvector::const_iterator i=add2begin; i!=add2end; ++i) { - add_indices = get_all_dummy_indices_safely(i->rest); - add2_dummy_indices.insert(add2_dummy_indices.end(), add_indices.begin(), add_indices.end()); - } + if (!skip_idx_rename) { + for (epvector::const_iterator i=add1begin; i!=add1end; ++i) { + add_indices = get_all_dummy_indices_safely(i->rest); + add1_dummy_indices.insert(add1_dummy_indices.end(), add_indices.begin(), add_indices.end()); + } + for (epvector::const_iterator i=add2begin; i!=add2end; ++i) { + add_indices = get_all_dummy_indices_safely(i->rest); + add2_dummy_indices.insert(add2_dummy_indices.end(), add_indices.begin(), add_indices.end()); + } - sort(add1_dummy_indices.begin(), add1_dummy_indices.end(), ex_is_less()); - sort(add2_dummy_indices.begin(), add2_dummy_indices.end(), ex_is_less()); - lst dummy_subs = rename_dummy_indices_uniquely(add1_dummy_indices, add2_dummy_indices); + sort(add1_dummy_indices.begin(), add1_dummy_indices.end(), ex_is_less()); + sort(add2_dummy_indices.begin(), add2_dummy_indices.end(), ex_is_less()); + dummy_subs = rename_dummy_indices_uniquely(add1_dummy_indices, add2_dummy_indices); + } // Multiply explicitly all non-numeric terms of add1 and add2: for (epvector::const_iterator i2=add2begin; i2!=add2end; ++i2) { // We really have to combine terms here in order to compactify // the result. Otherwise it would become waayy tooo bigg. - numeric oc; - distrseq.clear(); - ex i2_new = (dummy_subs.op(0).nops()>0? - i2->rest.subs((lst)dummy_subs.op(0), (lst)dummy_subs.op(1), subs_options::no_pattern) : i2->rest); + numeric oc(*_num0_p); + epvector distrseq2; + distrseq2.reserve(add1.seq.size()); + const ex i2_new = (skip_idx_rename || (dummy_subs.op(0).nops() == 0) ? + i2->rest : + i2->rest.subs(ex_to<lst>(dummy_subs.op(0)), + ex_to<lst>(dummy_subs.op(1)), subs_options::no_pattern)); for (epvector::const_iterator i1=add1begin; i1!=add1end; ++i1) { // Don't push_back expairs which might have a rest that evaluates to a numeric, // since that would violate an invariant of expairseq: @@ -1076,13 +1083,12 @@ ex mul::expand(unsigned options) const if (is_exactly_a<numeric>(rest)) { oc += ex_to<numeric>(rest).mul(ex_to<numeric>(i1->coeff).mul(ex_to<numeric>(i2->coeff))); } else { - distrseq.push_back(expair(rest, ex_to<numeric>(i1->coeff).mul_dyn(ex_to<numeric>(i2->coeff)))); + distrseq2.push_back(expair(rest, ex_to<numeric>(i1->coeff).mul_dyn(ex_to<numeric>(i2->coeff)))); } } - tmp_accu += (new add(distrseq, oc))->setflag(status_flags::dynallocated); - } + tmp_accu += (new add(distrseq2, oc))->setflag(status_flags::dynallocated); + } last_expanded = tmp_accu; - } else { if (!last_expanded.is_equal(_ex1)) non_adds.push_back(split_ex_to_pair(last_expanded)); @@ -1100,12 +1106,18 @@ ex mul::expand(unsigned options) const size_t n = last_expanded.nops(); exvector distrseq; distrseq.reserve(n); - exvector va = get_all_dummy_indices_safely(mul(non_adds)); - sort(va.begin(), va.end(), ex_is_less()); + exvector va; + if (! skip_idx_rename) { + va = get_all_dummy_indices_safely(mul(non_adds)); + sort(va.begin(), va.end(), ex_is_less()); + } for (size_t i=0; i<n; ++i) { epvector factors = non_adds; - factors.push_back(split_ex_to_pair(rename_dummy_indices_uniquely(va, last_expanded.op(i)))); + if (skip_idx_rename) + factors.push_back(split_ex_to_pair(last_expanded.op(i))); + else + factors.push_back(split_ex_to_pair(rename_dummy_indices_uniquely(va, last_expanded.op(i)))); ex term = (new mul(factors, overall_coeff))->setflag(status_flags::dynallocated); if (can_be_further_expanded(term)) { distrseq.push_back(term.expand()); -- 1.5.2.4 Best regards, Alexei -- All science is either physics or stamp collecting.
Apparently, in ~ 30% of calls to mul::expand the expression is monomial. Expanding monomials should be done as fast as possible. --- ginac/mul.cpp | 10 ++++++++++ ginac/power.cpp | 3 +++ 2 files changed, 13 insertions(+), 0 deletions(-) diff --git a/ginac/mul.cpp b/ginac/mul.cpp index db8b9f1..298715b 100644 --- a/ginac/mul.cpp +++ b/ginac/mul.cpp @@ -34,6 +34,7 @@ #include "lst.h" #include "archive.h" #include "utils.h" +#include "symbol.h" #include "compiler.h" namespace GiNaC { @@ -992,6 +993,15 @@ bool mul::can_be_further_expanded(const ex & e) ex mul::expand(unsigned options) const { + { + // trivial case: expanding the monomial (~ 30% of all calls) + epvector::const_iterator i = seq.begin(), seq_end = seq.end(); + while ((i != seq.end()) && is_a<symbol>(i->rest) && i->coeff.info(info_flags::integer)) + ++i; + if (i == seq_end) + return (new mul(*this))->setflag(status_flags::dynallocated | status_flags::expanded); + } + const bool skip_idx_rename = ! info(info_flags::has_indices); // First, expand the children std::auto_ptr<epvector> expanded_seqp = expandchildren(options); diff --git a/ginac/power.cpp b/ginac/power.cpp index fe7132a..cc003bc 100644 --- a/ginac/power.cpp +++ b/ginac/power.cpp @@ -771,6 +771,9 @@ tinfo_t power::return_type_tinfo() const ex power::expand(unsigned options) const { + if (is_a<symbol>(basis) && exponent.info(info_flags::integer)) + return (new power(*this))->setflag(status_flags::dynallocated | status_flags::expanded); + if (options == 0 && (flags & status_flags::expanded)) return *this; -- 1.5.2.4 Best regards, Alexei -- All science is either physics or stamp collecting.
This helps mul::expand() and friends to recognize objects which have no indices at all. --- ginac/flags.h | 3 ++- ginac/mul.cpp | 8 +++++++- ginac/power.cpp | 7 ++++++- 3 files changed, 15 insertions(+), 3 deletions(-) diff --git a/ginac/flags.h b/ginac/flags.h index 372b899..39a2113 100644 --- a/ginac/flags.h +++ b/ginac/flags.h @@ -30,7 +30,8 @@ class expand_options { public: enum { expand_indexed = 0x0001, ///< expands (a+b).i to a.i+b.i - expand_function_args = 0x0002 ///< expands the arguments of functions + expand_function_args = 0x0002, ///< expands the arguments of functions + expand_rename_idx = 0x0004 ///< used internally by mul::expand() }; }; diff --git a/ginac/mul.cpp b/ginac/mul.cpp index 298715b..4e95a63 100644 --- a/ginac/mul.cpp +++ b/ginac/mul.cpp @@ -1002,7 +1002,13 @@ ex mul::expand(unsigned options) const return (new mul(*this))->setflag(status_flags::dynallocated | status_flags::expanded); } - const bool skip_idx_rename = ! info(info_flags::has_indices); + // do not rename indices if the object has no indices at all + if ((!(options & expand_options::expand_rename_idx)) && + this->info(info_flags::has_indices)) + options |= expand_options::expand_rename_idx; + + const bool skip_idx_rename = !(options & expand_options::expand_rename_idx); + // First, expand the children std::auto_ptr<epvector> expanded_seqp = expandchildren(options); const epvector & expanded_seq = (expanded_seqp.get() ? *expanded_seqp : seq); diff --git a/ginac/power.cpp b/ginac/power.cpp index cc003bc..7d909e8 100644 --- a/ginac/power.cpp +++ b/ginac/power.cpp @@ -1009,8 +1009,13 @@ ex power::expand_mul(const mul & m, const numeric & n, unsigned options, bool fr return _ex1; } + // do not bother to rename indices if there are no any. + if ((!(options & expand_options::expand_rename_idx)) + && m.info(info_flags::has_indices)) + options |= expand_options::expand_rename_idx; // Leave it to multiplication since dummy indices have to be renamed - if (get_all_dummy_indices(m).size() > 0 && n.is_positive()) { + if ((options & expand_options::expand_rename_idx) && + (get_all_dummy_indices(m).size() > 0) && n.is_positive()) { ex result = m; exvector va = get_all_dummy_indices(m); sort(va.begin(), va.end(), ex_is_less()); -- 1.5.2.4 Best regards, Alexei -- All science is either physics or stamp collecting.
Do not bother to rename indices if object has no indices at all. --- ginac/expairseq.cpp | 15 ++++++++++++--- 1 files changed, 12 insertions(+), 3 deletions(-) diff --git a/ginac/expairseq.cpp b/ginac/expairseq.cpp index 79c741a..788f783 100644 --- a/ginac/expairseq.cpp +++ b/ginac/expairseq.cpp @@ -812,8 +812,8 @@ void expairseq::construct_from_2_ex(const ex &lh, const ex &rh) construct_from_2_ex_via_exvector(lh,rh); } else { #endif // EXPAIRSEQ_USE_HASHTAB - if(is_a<mul>(lh)) - { + if (is_a<mul>(lh) && lh.info(info_flags::has_indices) && + rh.info(info_flags::has_indices)) { ex newrh=rename_dummy_indices_uniquely(lh, rh); construct_from_2_expairseq(ex_to<expairseq>(lh), ex_to<expairseq>(newrh)); @@ -1057,6 +1057,7 @@ void expairseq::make_flat(const exvector &v) // and their cumulative number of operands int nexpairseqs = 0; int noperands = 0; + bool do_idx_rename = false; cit = v.begin(); while (cit!=v.end()) { @@ -1064,6 +1065,9 @@ void expairseq::make_flat(const exvector &v) ++nexpairseqs; noperands += ex_to<expairseq>(*cit).seq.size(); } + if (is_a<mul>(*this) && (!do_idx_rename) && + cit->info(info_flags::has_indices)) + do_idx_rename = true; ++cit; } @@ -1071,7 +1075,7 @@ void expairseq::make_flat(const exvector &v) seq.reserve(v.size()+noperands-nexpairseqs); // copy elements and split off numerical part - make_flat_inserter mf(v, this->tinfo() == &mul::tinfo_static); + make_flat_inserter mf(v, do_idx_rename); cit = v.begin(); while (cit!=v.end()) { if (ex_to<basic>(*cit).tinfo()==this->tinfo()) { @@ -1105,6 +1109,7 @@ void expairseq::make_flat(const epvector &v, bool do_index_renaming) // and their cumulative number of operands int nexpairseqs = 0; int noperands = 0; + bool really_need_rename_inds = false; cit = v.begin(); while (cit!=v.end()) { @@ -1112,8 +1117,12 @@ void expairseq::make_flat(const epvector &v, bool do_index_renaming) ++nexpairseqs; noperands += ex_to<expairseq>(cit->rest).seq.size(); } + if ((!really_need_rename_inds) && is_a<mul>(*this) && + cit->rest.info(info_flags::has_indices)) + really_need_rename_inds = true; ++cit; } + do_index_renaming = do_index_renaming && really_need_rename_inds; // reserve seq and coeffseq which will hold all operands seq.reserve(v.size()+noperands-nexpairseqs); -- 1.5.2.4 Best regards, Alexei -- All science is either physics or stamp collecting.
This is the similar set of patches for GiNaC 1.3 --- ginac/expairseq.cpp | 23 +++++++++++++++++++++-- ginac/flags.h | 4 +++- ginac/idx.cpp | 7 +++++-- ginac/power.cpp | 15 +++++++++++++++ ginac/symbol.cpp | 30 ++++++++++++++++-------------- 5 files changed, 60 insertions(+), 19 deletions(-) diff --git a/ginac/expairseq.cpp b/ginac/expairseq.cpp index be8790f..97cf961 100644 --- a/ginac/expairseq.cpp +++ b/ginac/expairseq.cpp @@ -24,6 +24,7 @@ #include <algorithm> #include <string> #include <stdexcept> +#include <iterator> #include "expairseq.h" #include "lst.h" @@ -261,8 +262,26 @@ void expairseq::do_print_tree(const print_tree & c, unsigned level) const bool expairseq::info(unsigned inf) const { - if (inf == info_flags::expanded) - return (flags & status_flags::expanded); + switch(inf) { + case info_flags::expanded: + return (flags & status_flags::expanded); + case info_flags::has_indices: { + if (flags & status_flags::has_indices) + return true; + else if (flags & status_flags::has_no_indices) + return false; + for (epvector::const_iterator i = seq.begin(); i != seq.end(); ++i) { + if (i->rest.info(info_flags::has_indices)) { + this->setflag(status_flags::has_indices); + this->clearflag(status_flags::has_no_indices); + return true; + } + } + this->clearflag(status_flags::has_indices); + this->setflag(status_flags::has_no_indices); + return false; + } + } return inherited::info(inf); } diff --git a/ginac/flags.h b/ginac/flags.h index a42f869..33622f9 100644 --- a/ginac/flags.h +++ b/ginac/flags.h @@ -179,7 +179,9 @@ public: evaluated = 0x0002, ///< .eval() has already done its job expanded = 0x0004, ///< .expand(0) has already done its job (other expand() options ignore this flag) hash_calculated = 0x0008, ///< .calchash() has already done its job - not_shareable = 0x0010 ///< don't share instances of this object between different expressions unless explicitly asked to (used by ex::compare()) + not_shareable = 0x0010, ///< don't share instances of this object between different expressions unless explicitly asked to (used by ex::compare()) + has_indices = 0x0020, + has_no_indices = 0x0040 // ! (has_indices || has_no_indices) means "don't know" }; }; diff --git a/ginac/idx.cpp b/ginac/idx.cpp index 3d363e0..601a77c 100644 --- a/ginac/idx.cpp +++ b/ginac/idx.cpp @@ -232,8 +232,11 @@ void spinidx::do_print_tree(const print_tree & c, unsigned level) const bool idx::info(unsigned inf) const { - if (inf == info_flags::idx) - return true; + switch(inf) { + case info_flags::idx: + case info_flags::has_indices: + return true; + } return inherited::info(inf); } diff --git a/ginac/power.cpp b/ginac/power.cpp index f243147..676b862 100644 --- a/ginac/power.cpp +++ b/ginac/power.cpp @@ -240,6 +240,21 @@ bool power::info(unsigned inf) const basis.info(inf); case info_flags::expanded: return (flags & status_flags::expanded); + case info_flags::has_indices: { + if (flags & status_flags::has_indices) + return true; + else if (flags & status_flags::has_no_indices) + return false; + else if (basis.info(info_flags::has_indices)) { + setflag(status_flags::has_indices); + clearflag(status_flags::has_no_indices); + return true; + } else { + clearflag(status_flags::has_indices); + setflag(status_flags::has_no_indices); + return false; + } + } } return inherited::info(inf); } diff --git a/ginac/symbol.cpp b/ginac/symbol.cpp index 2edd714..0e51f0e 100644 --- a/ginac/symbol.cpp +++ b/ginac/symbol.cpp @@ -190,20 +190,22 @@ void symbol::do_print_python_repr(const print_python_repr & c, unsigned level) c bool symbol::info(unsigned inf) const { - if (inf == info_flags::symbol) - return true; - if (inf == info_flags::polynomial || - inf == info_flags::integer_polynomial || - inf == info_flags::cinteger_polynomial || - inf == info_flags::rational_polynomial || - inf == info_flags::crational_polynomial || - inf == info_flags::rational_function || - inf == info_flags::expanded) - return true; - if (inf == info_flags::real) - return domain == domain::real; - else - return inherited::info(inf); + switch (inf) { + case info_flags::symbol: + case info_flags::polynomial: + case info_flags::integer_polynomial: + case info_flags::cinteger_polynomial: + case info_flags::rational_polynomial: + case info_flags::crational_polynomial: + case info_flags::rational_function: + case info_flags::expanded: + return true; + case info_flags::has_indices: + return false; + case info_flags::real: + return domain == domain::real; + } + return inherited::info(inf); } ex symbol::eval(int level) const -- 1.5.2.4 Best regards, Alexei -- All science is either physics or stamp collecting.
--- ginac/mul.cpp | 8 ++++++-- 1 files changed, 6 insertions(+), 2 deletions(-) diff --git a/ginac/mul.cpp b/ginac/mul.cpp index 3c87d5a..5176218 100644 --- a/ginac/mul.cpp +++ b/ginac/mul.cpp @@ -940,6 +940,7 @@ bool mul::can_be_further_expanded(const ex & e) ex mul::expand(unsigned options) const { + const bool skip_idx_rename = ! info(info_flags::has_indices); // First, expand the children std::auto_ptr<epvector> expanded_seqp = expandchildren(options); const epvector & expanded_seq = (expanded_seqp.get() ? *expanded_seqp : seq); @@ -1003,7 +1004,7 @@ ex mul::expand(unsigned options) const for (epvector::const_iterator i2=add2begin; i2!=add2end; ++i2) { // Don't push_back expairs which might have a rest that evaluates to a numeric, // since that would violate an invariant of expairseq: - const ex rest = (new mul(i1->rest, rename_dummy_indices_uniquely(i1->rest, i2->rest)))->setflag(status_flags::dynallocated); + const ex rest = (new mul(i1->rest, skip_idx_rename ? i2->rest : rename_dummy_indices_uniquely(i1->rest, i2->rest)))->setflag(status_flags::dynallocated); if (is_exactly_a<numeric>(rest)) { oc += ex_to<numeric>(rest).mul(ex_to<numeric>(i1->coeff).mul(ex_to<numeric>(i2->coeff))); } else { @@ -1034,7 +1035,10 @@ ex mul::expand(unsigned options) const for (size_t i=0; i<n; ++i) { epvector factors = non_adds; - factors.push_back(split_ex_to_pair(rename_dummy_indices_uniquely(mul(non_adds), last_expanded.op(i)))); + if (skip_idx_rename) + factors.push_back(split_ex_to_pair(last_expanded.op(i))); + else + factors.push_back(split_ex_to_pair(rename_dummy_indices_uniquely(mul(non_adds), last_expanded.op(i)))); ex term = (new mul(factors, overall_coeff))->setflag(status_flags::dynallocated); if (can_be_further_expanded(term)) { distrseq.push_back(term.expand()); -- 1.5.2.4 Best regards, Alexei -- All science is either physics or stamp collecting.
--- ginac/power.cpp | 2 +- 1 files changed, 1 insertions(+), 1 deletions(-) diff --git a/ginac/power.cpp b/ginac/power.cpp index 676b862..bbcb453 100644 --- a/ginac/power.cpp +++ b/ginac/power.cpp @@ -916,7 +916,7 @@ ex power::expand_mul(const mul & m, const numeric & n, unsigned options, bool fr } // Leave it to multiplication since dummy indices have to be renamed - if (get_all_dummy_indices(m).size() > 0 && n.is_positive()) { + if (m.info(info_flags::has_indices) && (get_all_dummy_indices(m).size() > 0) && n.is_positive()) { ex result = m; for (int i=1; i < n.to_int(); i++) result *= rename_dummy_indices_uniquely(m,m); -- 1.5.2.4 Best regards, Alexei -- All science is either physics or stamp collecting.
--- ginac/mul.cpp | 12 ++++++++++++ ginac/power.cpp | 3 +++ 2 files changed, 15 insertions(+), 0 deletions(-) diff --git a/ginac/mul.cpp b/ginac/mul.cpp index 5176218..d6ff7cb 100644 --- a/ginac/mul.cpp +++ b/ginac/mul.cpp @@ -34,6 +34,7 @@ #include "lst.h" #include "archive.h" #include "utils.h" +#include "symbol.h" #include "compiler.h" namespace GiNaC { @@ -940,7 +941,18 @@ bool mul::can_be_further_expanded(const ex & e) ex mul::expand(unsigned options) const { + { + // trivial case: expanding the monomial (~ 30% of all calls) + epvector::const_iterator i = seq.begin(), seq_end = seq.end(); + while ((i != seq.end()) && is_a<symbol>(i->rest) && i->coeff.info(info_flags::integer)) + ++i; + if (i == seq_end) + return (new mul(*this))->setflag(status_flags::dynallocated | status_flags::expanded); + } + const bool skip_idx_rename = ! info(info_flags::has_indices); + if (skip_idx_rename) + ++(mul_expand_stats::n_indexless); // First, expand the children std::auto_ptr<epvector> expanded_seqp = expandchildren(options); const epvector & expanded_seq = (expanded_seqp.get() ? *expanded_seqp : seq); diff --git a/ginac/power.cpp b/ginac/power.cpp index bbcb453..010af6c 100644 --- a/ginac/power.cpp +++ b/ginac/power.cpp @@ -680,6 +680,9 @@ unsigned power::return_type_tinfo() const ex power::expand(unsigned options) const { + if (is_a<symbol>(basis) && exponent.info(info_flags::integer)) + return (new power(*this))->setflag(status_flags::dynallocated | status_flags::expanded); + if (options == 0 && (flags & status_flags::expanded)) return *this; -- 1.5.2.4 Best regards, Alexei -- All science is either physics or stamp collecting.
Helps to avoid useless (and expansive) rename_dummy_indices_uniquely(). --- ginac/flags.h | 3 ++- ginac/mul.cpp | 10 +++++++--- ginac/power.cpp | 7 ++++++- 3 files changed, 15 insertions(+), 5 deletions(-) diff --git a/ginac/flags.h b/ginac/flags.h index 33622f9..478da32 100644 --- a/ginac/flags.h +++ b/ginac/flags.h @@ -30,7 +30,8 @@ class expand_options { public: enum { expand_indexed = 0x0001, ///< expands (a+b).i to a.i+b.i - expand_function_args = 0x0002 ///< expands the arguments of functions + expand_function_args = 0x0002, ///< expands the arguments of functions + expand_rename_idx = 0x0004 ///< used internally by mul::expand() }; }; diff --git a/ginac/mul.cpp b/ginac/mul.cpp index d6ff7cb..b7e0e82 100644 --- a/ginac/mul.cpp +++ b/ginac/mul.cpp @@ -950,9 +950,13 @@ ex mul::expand(unsigned options) const return (new mul(*this))->setflag(status_flags::dynallocated | status_flags::expanded); } - const bool skip_idx_rename = ! info(info_flags::has_indices); - if (skip_idx_rename) - ++(mul_expand_stats::n_indexless); + // do not rename indices if the object has no indices at all + if ((!(options & expand_options::expand_rename_idx)) && + this->info(info_flags::has_indices)) + options |= expand_options::expand_rename_idx; + + const bool skip_idx_rename = !(options & expand_options::expand_rename_idx); + // First, expand the children std::auto_ptr<epvector> expanded_seqp = expandchildren(options); const epvector & expanded_seq = (expanded_seqp.get() ? *expanded_seqp : seq); diff --git a/ginac/power.cpp b/ginac/power.cpp index 010af6c..42e14c1 100644 --- a/ginac/power.cpp +++ b/ginac/power.cpp @@ -918,8 +918,13 @@ ex power::expand_mul(const mul & m, const numeric & n, unsigned options, bool fr return _ex1; } + // do not bother to rename indices if there are no any. + if ((!(options & expand_options::expand_rename_idx)) + && m.info(info_flags::has_indices)) + options |= expand_options::expand_rename_idx; // Leave it to multiplication since dummy indices have to be renamed - if (m.info(info_flags::has_indices) && (get_all_dummy_indices(m).size() > 0) && n.is_positive()) { + if ((!(options & expand_options::expand_rename_idx)) && + (get_all_dummy_indices(m).size() > 0) && n.is_positive()) { ex result = m; for (int i=1; i < n.to_int(); i++) result *= rename_dummy_indices_uniquely(m,m); -- 1.5.2.4 -- All science is either physics or stamp collecting.
Alexei Sheplyakov wrote:
This is the similar set of patches for GiNaC 1.3
A simple proposal: let's not care about release 1.3 so much. It is a waste of time now that 1.4 is out. If there's some really serious issue with 1.3 we can fix it for the 1.3 branch, but only as long as it is still in use by people out there. But making an effort to backport as much as possible appears pointless to me. And as soon as there are radical proposals that break the ABI it would make more sense to branch 1.4. Cheers -richy. -- Richard B. Kreckel <http://www.ginac.de/~kreckel/>
Hello! On Mon, Sep 17, 2007 at 09:30:04PM +0200, Richard B. Kreckel wrote:
A simple proposal: let's not care about release 1.3 so much. It is a waste of time now that 1.4 is out.
Frankly, I'm not going to switch to 1.4 anytime soon, so I do care about 1.3. In general, I don't like to upgrade unless there is a compelling reason to do so (e.g. new version is 2x faster, or has some feature I need very much, etc).
If there's some really serious issue with 1.3 we can fix it for the 1.3 branch,
I think this particular regression is pretty serious.
but only as long as it is still in use by people out there.
There is at least 1 user :) People don't like to upgrade (especially if their code breaks, gets substantially slower, etc), so I think there will be more.
But making an effort to backport as much as possible appears pointless to me.
Actually, [almost] all patches for 1.4 (and HEAD) I've submitted were forward ported (in most cases porting boils down to git-cherry-pick). Best regards, Alexei -- All science is either physics or stamp collecting.
participants (3)
-
Alexei Sheplyakov
-
Richard B. Kreckel
-
varg＠theor.jinr.ru