X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=blobdiff_plain;f=ginac%2Fmul.cpp;h=27add930cf56863ae7d635d3031e4dc2b1d7ea9f;hp=fe63a08861f594b3ca79f6e421879a6d235c137a;hb=f320e27e9cfe287168c879af5991babffaa1e9c8;hpb=96af2609db413ddf70371e0423e6c07ecb5ee813 diff --git a/ginac/mul.cpp b/ginac/mul.cpp index fe63a088..27add930 100644 --- a/ginac/mul.cpp +++ b/ginac/mul.cpp @@ -3,7 +3,7 @@ * Implementation of GiNaC's products of expressions. */ /* - * GiNaC Copyright (C) 1999-2004 Johannes Gutenberg University Mainz, Germany + * GiNaC Copyright (C) 1999-2011 Johannes Gutenberg University Mainz, Germany * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by @@ -17,22 +17,25 @@ * * You should have received a copy of the GNU General Public License * along with this program; if not, write to the Free Software - * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA */ -#include -#include -#include -#include - #include "mul.h" #include "add.h" #include "power.h" #include "operators.h" #include "matrix.h" +#include "indexed.h" #include "lst.h" #include "archive.h" #include "utils.h" +#include "symbol.h" +#include "compiler.h" + +#include +#include +#include +#include namespace GiNaC { @@ -50,7 +53,6 @@ GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(mul, expairseq, mul::mul() { - tinfo_key = TINFO_mul; } ////////// @@ -61,7 +63,6 @@ mul::mul() mul::mul(const ex & lh, const ex & rh) { - tinfo_key = TINFO_mul; overall_coeff = _ex1; construct_from_2_ex(lh,rh); GINAC_ASSERT(is_canonical()); @@ -69,7 +70,6 @@ mul::mul(const ex & lh, const ex & rh) mul::mul(const exvector & v) { - tinfo_key = TINFO_mul; overall_coeff = _ex1; construct_from_exvector(v); GINAC_ASSERT(is_canonical()); @@ -77,32 +77,28 @@ mul::mul(const exvector & v) mul::mul(const epvector & v) { - tinfo_key = TINFO_mul; overall_coeff = _ex1; construct_from_epvector(v); GINAC_ASSERT(is_canonical()); } -mul::mul(const epvector & v, const ex & oc) +mul::mul(const epvector & v, const ex & oc, bool do_index_renaming) { - tinfo_key = TINFO_mul; overall_coeff = oc; - construct_from_epvector(v); + construct_from_epvector(v, do_index_renaming); GINAC_ASSERT(is_canonical()); } -mul::mul(std::auto_ptr vp, const ex & oc) +mul::mul(std::auto_ptr vp, const ex & oc, bool do_index_renaming) { - tinfo_key = TINFO_mul; - GINAC_ASSERT(vp!=0); + GINAC_ASSERT(vp.get()!=0); overall_coeff = oc; - construct_from_epvector(*vp); + construct_from_epvector(*vp, do_index_renaming); GINAC_ASSERT(is_canonical()); } mul::mul(const ex & lh, const ex & mh, const ex & rh) { - tinfo_key = TINFO_mul; exvector factors; factors.reserve(3); factors.push_back(lh); @@ -117,8 +113,6 @@ mul::mul(const ex & lh, const ex & mh, const ex & rh) // archiving ////////// -DEFAULT_ARCHIVING(mul) - ////////// // functions overriding virtual functions from base classes ////////// @@ -128,8 +122,8 @@ void mul::print_overall_coeff(const print_context & c, const char *mul_sym) cons const numeric &coeff = ex_to(overall_coeff); if (coeff.csgn() == -1) c.s << '-'; - if (!coeff.is_equal(_num1) && - !coeff.is_equal(_num_1)) { + if (!coeff.is_equal(*_num1_p) && + !coeff.is_equal(*_num_1_p)) { if (coeff.is_rational()) { if (coeff.is_negative()) (-coeff).print(c); @@ -217,8 +211,12 @@ void mul::do_print_csrc(const print_csrc & c, unsigned level) const c.s << "("; if (!overall_coeff.is_equal(_ex1)) { - overall_coeff.print(c, precedence()); - c.s << "*"; + if (overall_coeff.is_equal(_ex_1)) + c.s << "-"; + else { + overall_coeff.print(c, precedence()); + c.s << "*"; + } } // Print arguments, separated by "*" or "/" @@ -280,6 +278,16 @@ bool mul::info(unsigned inf) const case info_flags::integer_polynomial: case info_flags::cinteger_polynomial: case info_flags::rational_polynomial: + case info_flags::real: + case info_flags::rational: + case info_flags::integer: + case info_flags::crational: + case info_flags::cinteger: + case info_flags::positive: + case info_flags::nonnegative: + case info_flags::posint: + case info_flags::nonnegint: + case info_flags::even: case info_flags::crational_polynomial: case info_flags::rational_function: { epvector::const_iterator i = seq.begin(), end = seq.end(); @@ -288,6 +296,8 @@ bool mul::info(unsigned inf) const return false; ++i; } + if (overall_coeff.is_equal(*_num1_p) && inf == info_flags::even) + return true; return overall_coeff.info(inf); } case info_flags::algebraic: { @@ -299,10 +309,55 @@ bool mul::info(unsigned inf) const } return false; } + case info_flags::negative: { + bool neg = false; + epvector::const_iterator i = seq.begin(), end = seq.end(); + while (i != end) { + const ex& factor = recombine_pair_to_ex(*i++); + if (factor.info(info_flags::positive)) + continue; + else if (factor.info(info_flags::negative)) + neg = !neg; + else + return false; + } + if (overall_coeff.info(info_flags::negative)) + neg = !neg; + return neg; + } + case info_flags::negint: { + bool neg = false; + epvector::const_iterator i = seq.begin(), end = seq.end(); + while (i != end) { + const ex& factor = recombine_pair_to_ex(*i++); + if (factor.info(info_flags::posint)) + continue; + else if (factor.info(info_flags::negint)) + neg = !neg; + else + return false; + } + if (overall_coeff.info(info_flags::negint)) + neg = !neg; + else if (!overall_coeff.info(info_flags::posint)) + return false; + return neg; + } } return inherited::info(inf); } +bool mul::is_polynomial(const ex & var) const +{ + for (epvector::const_iterator i=seq.begin(); i!=seq.end(); ++i) { + if (!i->rest.is_polynomial(var) || + (i->rest.has(var) && !i->coeff.info(info_flags::integer))) { + return false; + } + } + return true; +} + int mul::degree(const ex & s) const { // Sum up degrees of factors @@ -310,7 +365,11 @@ int mul::degree(const ex & s) const epvector::const_iterator i = seq.begin(), end = seq.end(); while (i != end) { if (ex_to(i->coeff).is_integer()) - deg_sum += i->rest.degree(s) * ex_to(i->coeff).to_int(); + deg_sum += recombine_pair_to_ex(*i).degree(s); + else { + if (i->rest.has(s)) + throw std::runtime_error("mul::degree() undefined degree because of non-integer exponent"); + } ++i; } return deg_sum; @@ -323,7 +382,11 @@ int mul::ldegree(const ex & s) const epvector::const_iterator i = seq.begin(), end = seq.end(); while (i != end) { if (ex_to(i->coeff).is_integer()) - deg_sum += i->rest.ldegree(s) * ex_to(i->coeff).to_int(); + deg_sum += recombine_pair_to_ex(*i).ldegree(s); + else { + if (i->rest.has(s)) + throw std::runtime_error("mul::ldegree() undefined degree because of non-integer exponent"); + } ++i; } return deg_sum; @@ -409,7 +472,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; @@ -421,7 +484,7 @@ ex mul::eval(int level) const return recombine_pair_to_ex(*(seq.begin())); } else if ((seq_size==1) && is_exactly_a((*seq.begin()).rest) && - ex_to((*seq.begin()).coeff).is_equal(_num1)) { + ex_to((*seq.begin()).coeff).is_equal(*_num1_p)) { // *(+(x,y,...);c) -> +(*(x,c),*(y,c),...) (c numeric(), no powers of +()) const add & addref = ex_to((*seq.begin()).rest); std::auto_ptr distrseq(new epvector); @@ -433,9 +496,81 @@ ex mul::eval(int level) const } return (new add(distrseq, ex_to(addref.overall_coeff). - mul_dyn(ex_to(overall_coeff)))) - ->setflag(status_flags::dynallocated | status_flags::evaluated); + mul_dyn(ex_to(overall_coeff))) + )->setflag(status_flags::dynallocated | status_flags::evaluated); + } 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 s(new epvector); + numeric oc = *_num1_p; + bool something_changed = false; + while (i!=last) { + if (likely(! (is_a(i->rest) && i->coeff.is_equal(_ex1)))) { + // power::eval has such a rule, no need to handle powers here + ++i; + continue; + } + + // XXX: What is the best way to check if the polynomial is a primitive? + numeric c = i->rest.integer_content(); + const numeric lead_coeff = + ex_to(ex_to(i->rest).seq.begin()->coeff).div(c); + const bool canonicalizable = lead_coeff.is_integer(); + + // XXX: The main variable is chosen in a random way, so this code + // does NOT transform the term into the canonical form (thus, in some + // 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 (likely((c == *_num1_p) && ((! canonicalizable) || unit_normal))) { + ++i; + continue; + } + + if (! something_changed) { + s->reserve(seq_size); + something_changed = true; + } + + while ((j!=i) && (j!=last)) { + s->push_back(*j); + ++j; + } + + if (! unit_normal) + c = c.mul(*_num_1_p); + + oc = oc.mul(c); + + // divide add by the number in place to save at least 2 .eval() calls + const add& addref = ex_to(i->rest); + add* primitive = new add(addref); + primitive->setflag(status_flags::dynallocated); + primitive->clearflag(status_flags::hash_calculated); + primitive->overall_coeff = ex_to(primitive->overall_coeff).div_dyn(c); + for (epvector::iterator ai = primitive->seq.begin(); + ai != primitive->seq.end(); ++ai) + ai->coeff = ex_to(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(overall_coeff).mul_dyn(oc)) + )->setflag(status_flags::dynallocated); + } } + return this->hold(); } @@ -460,6 +595,41 @@ ex mul::evalf(int level) const return mul(s, overall_coeff.evalf(level)); } +void mul::find_real_imag(ex & rp, ex & ip) const +{ + rp = overall_coeff.real_part(); + ip = overall_coeff.imag_part(); + for (epvector::const_iterator i=seq.begin(); i!=seq.end(); ++i) { + ex factor = recombine_pair_to_ex(*i); + ex new_rp = factor.real_part(); + ex new_ip = factor.imag_part(); + if(new_ip.is_zero()) { + rp *= new_rp; + ip *= new_rp; + } else { + ex temp = rp*new_rp - ip*new_ip; + ip = ip*new_rp + rp*new_ip; + rp = temp; + } + } + rp = rp.expand(); + ip = ip.expand(); +} + +ex mul::real_part() const +{ + ex rp, ip; + find_real_imag(rp, ip); + return rp; +} + +ex mul::imag_part() const +{ + ex rp, ip; + find_real_imag(rp, ip); + return ip; +} + ex mul::evalm() const { // numeric*matrix @@ -515,7 +685,7 @@ ex mul::eval_ncmul(const exvector & v) const return inherited::eval_ncmul(v); } -bool tryfactsubs(const ex & origfactor, const ex & patternfactor, int & nummatches, lst & repls) +bool tryfactsubs(const ex & origfactor, const ex & patternfactor, int & nummatches, exmap& repls) { ex origbase; int origexponent; @@ -547,7 +717,7 @@ bool tryfactsubs(const ex & origfactor, const ex & patternfactor, int & nummatch patternexpsign = 1; } - lst saverepls = repls; + exmap saverepls = repls; if (origexponent < patternexponent || origexpsign != patternexpsign || !origbase.match(patternbase,saverepls)) return false; repls = saverepls; @@ -558,61 +728,102 @@ bool tryfactsubs(const ex & origfactor, const ex & patternfactor, int & nummatch return true; } +/** Checks wheter e matches to the pattern pat and the (possibly to be updated) + * list of replacements repls. This matching is in the sense of algebraic + * substitutions. Matching starts with pat.op(factor) of the pattern because + * the factors before this one have already been matched. The (possibly + * updated) number of matches is in nummatches. subsed[i] is true for factors + * that already have been replaced by previous substitutions and matched[i] + * is true for factors that have been matched by the current match. + */ +bool algebraic_match_mul_with_mul(const mul &e, const ex &pat, exmap& repls, + int factor, int &nummatches, const std::vector &subsed, + std::vector &matched) +{ + GINAC_ASSERT(subsed.size() == e.nops()); + GINAC_ASSERT(matched.size() == e.nops()); + + if (factor == (int)pat.nops()) + return true; + + for (size_t i=0; i(pattern)) { + exmap repls; + int nummatches = std::numeric_limits::max(); + std::vector subsed(nops(), false); + std::vector matched(nops(), false); + if(algebraic_match_mul_with_mul(*this, pattern, repls, 0, nummatches, + subsed, matched)) + return true; + } + return basic::has(pattern, options); +} + ex mul::algebraic_subs_mul(const exmap & m, unsigned options) const { - std::vector subsed(seq.size(), false); - exvector subsresult(seq.size()); + std::vector subsed(nops(), false); + ex divide_by = 1; + ex multiply_by = 1; for (exmap::const_iterator it = m.begin(); it != m.end(); ++it) { if (is_exactly_a(it->first)) { - +retry1: int nummatches = std::numeric_limits::max(); - std::vector currsubsed(seq.size(), false); - bool succeed = true; - lst repls; - - for (size_t j=0; jfirst.nops(); j++) { - bool found=false; - for (size_t k=0; kfirst.op(j), nummatches, repls)) { - currsubsed[k] = true; - found = true; - break; - } - } - if (!found) { - succeed = false; - break; - } - } - if (!succeed) + std::vector currsubsed(nops(), false); + exmap repls; + + if(!algebraic_match_mul_with_mul(*this, it->first, repls, 0, nummatches, subsed, currsubsed)) continue; - bool foundfirstsubsedfactor = false; - for (size_t j=0; jsecond.subs(ex(repls), subs_options::no_pattern) / it->first.subs(ex(repls), subs_options::no_pattern), nummatches); - } + for (size_t j=0; jfirst.subs(repls, subs_options::no_pattern); + divide_by *= power(subsed_pattern, nummatches); + ex subsed_result + = it->second.subs(repls, subs_options::no_pattern); + multiply_by *= power(subsed_result, nummatches); + goto retry1; } else { - int nummatches = std::numeric_limits::max(); - lst repls; - for (size_t j=0; jnops(); j++) { - if (!subsed[j] && tryfactsubs(op(j), it->first, nummatches, repls)) { + int nummatches = std::numeric_limits::max(); + exmap repls; + if (!subsed[j] && tryfactsubs(op(j), it->first, nummatches, repls)){ subsed[j] = true; - subsresult[j] = op(j) * power(it->second.subs(ex(repls), subs_options::no_pattern) / it->first.subs(ex(repls), subs_options::no_pattern), nummatches); + ex subsed_pattern + = it->first.subs(repls, subs_options::no_pattern); + divide_by *= power(subsed_pattern, nummatches); + ex subsed_result + = it->second.subs(repls, subs_options::no_pattern); + multiply_by *= power(subsed_result, nummatches); } } } @@ -628,17 +839,41 @@ ex mul::algebraic_subs_mul(const exmap & m, unsigned options) const if (!subsfound) return subs_one_level(m, options | subs_options::algebraic); - exvector ev; ev.reserve(nops()); - for (size_t i=0; isetflag(status_flags::dynallocated); +ex mul::conjugate() const +{ + // The base class' method is wrong here because we have to be careful at + // branch cuts. power::conjugate takes care of that already, so use it. + epvector *newepv = 0; + for (epvector::const_iterator i=seq.begin(); i!=seq.end(); ++i) { + if (newepv) { + newepv->push_back(split_ex_to_pair(recombine_pair_to_ex(*i).conjugate())); + continue; + } + ex x = recombine_pair_to_ex(*i); + ex c = x.conjugate(); + if (c.is_equal(x)) { + continue; + } + newepv = new epvector; + newepv->reserve(seq.size()); + for (epvector::const_iterator j=seq.begin(); j!=i; ++j) { + newepv->push_back(*j); + } + newepv->push_back(split_ex_to_pair(c)); + } + ex x = overall_coeff.conjugate(); + if (!newepv && are_ex_trivially_equal(x, overall_coeff)) { + return *this; + } + ex result = thisexpairseq(newepv ? *newepv : seq, x); + delete newepv; + return result; } + // protected /** Implementation of ex::diff() for a product. It applies the product rule. @@ -672,7 +907,7 @@ int mul::compare_same_type(const basic & other) const unsigned mul::return_type() const { if (seq.empty()) { - // mul without factors: should not happen, but commutes + // mul without factors: should not happen, but commutates return return_types::commutative; } @@ -692,8 +927,8 @@ unsigned mul::return_type() const if ((rt == return_types::noncommutative) && (!all_commutative)) { // another nc element found, compare type_infos if (noncommutative_element->rest.return_type_tinfo() != i->rest.return_type_tinfo()) { - // diffent types -> mul is ncc - return return_types::noncommutative_composite; + // different types -> mul is ncc + return return_types::noncommutative_composite; } } ++i; @@ -702,10 +937,10 @@ unsigned mul::return_type() const return all_commutative ? return_types::commutative : return_types::noncommutative; } -unsigned mul::return_type_tinfo() const +return_type_t mul::return_type_tinfo() const { if (seq.empty()) - return tinfo_key; // mul without factors: should not happen + return make_return_type_t(); // mul without factors: should not happen // return type_info of first noncommutative element epvector::const_iterator i = seq.begin(), end = seq.end(); @@ -715,17 +950,17 @@ unsigned mul::return_type_tinfo() const ++i; } // no noncommutative element found, should not happen - return tinfo_key; + return make_return_type_t(); } -ex mul::thisexpairseq(const epvector & v, const ex & oc) const +ex mul::thisexpairseq(const epvector & v, const ex & oc, bool do_index_renaming) const { - return (new mul(v, oc))->setflag(status_flags::dynallocated); + return (new mul(v, oc, do_index_renaming))->setflag(status_flags::dynallocated); } -ex mul::thisexpairseq(std::auto_ptr vp, const ex & oc) const +ex mul::thisexpairseq(std::auto_ptr vp, const ex & oc, bool do_index_renaming) const { - return (new mul(vp, oc))->setflag(status_flags::dynallocated); + return (new mul(vp, oc, do_index_renaming))->setflag(status_flags::dynallocated); } expair mul::split_ex_to_pair(const ex & e) const @@ -766,7 +1001,7 @@ expair mul::combine_pair_with_coeff_to_pair(const expair & p, ex mul::recombine_pair_to_ex(const expair & p) const { - if (ex_to(p.coeff).is_equal(_num1)) + if (ex_to(p.coeff).is_equal(*_num1_p)) return p.rest; else return (new power(p.rest,p.coeff))->setflag(status_flags::dynallocated); @@ -821,7 +1056,7 @@ bool mul::can_make_flat(const expair & p) const // this assertion will probably fail somewhere // it would require a more careful make_flat, obeying the power laws // probably should return true only if p.coeff is integer - return ex_to(p.coeff).is_equal(_num1); + return ex_to(p.coeff).is_equal(*_num1_p); } bool mul::can_be_further_expanded(const ex & e) @@ -840,6 +1075,24 @@ 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(i->rest) && i->coeff.info(info_flags::integer)) + ++i; + if (i == seq_end) { + setflag(status_flags::expanded); + return *this; + } + } + + // 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 expanded_seqp = expandchildren(options); const epvector & expanded_seq = (expanded_seqp.get() ? *expanded_seqp : seq); @@ -848,7 +1101,6 @@ ex mul::expand(unsigned options) const // with the next one that is found while collecting the factors which are // not sums ex last_expanded = _ex1; - bool need_reexpand = false; epvector non_adds; non_adds.reserve(expanded_seq.size()); @@ -895,25 +1147,48 @@ ex mul::expand(unsigned options) const // Compute the new overall coefficient and put it together: 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; + + 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()); + 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 i1=add1begin; i1!=add1end; ++i1) { + 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(); - for (epvector::const_iterator i2=add2begin; i2!=add2end; ++i2) { + 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(dummy_subs.op(0)), + ex_to(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: - const ex rest = (new mul(i1->rest, i2->rest))->setflag(status_flags::dynallocated); - if (is_exactly_a(rest)) + const ex rest = (new mul(i1->rest, i2_new))->setflag(status_flags::dynallocated); + if (is_exactly_a(rest)) { oc += ex_to(rest).mul(ex_to(i1->coeff).mul(ex_to(i2->coeff))); - else - distrseq.push_back(expair(rest, ex_to(i1->coeff).mul_dyn(ex_to(i2->coeff)))); + } else { + distrseq2.push_back(expair(rest, ex_to(i1->coeff).mul_dyn(ex_to(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)); @@ -931,14 +1206,22 @@ ex mul::expand(unsigned options) const size_t n = last_expanded.nops(); exvector distrseq; distrseq.reserve(n); + 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; isetflag(status_flags::dynallocated); - if (can_be_further_expanded(term)) + if (can_be_further_expanded(term)) { distrseq.push_back(term.expand()); - else { + } else { if (options == 0) ex_to(term).setflag(status_flags::expanded); distrseq.push_back(term); @@ -1016,4 +1299,6 @@ std::auto_ptr mul::expandchildren(unsigned options) const return std::auto_ptr(0); // nothing has changed } +GINAC_BIND_UNARCHIVER(mul); + } // namespace GiNaC