X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=blobdiff_plain;f=ginac%2Fmul.cpp;h=16bbc5b229ac9f21c26bcadbfacb1a26c5b6790f;hp=9ed27c6c4bf3f4c7435674de9138bb4aa16d76fc;hb=f7884835d397de85e648d1957c058b7d4c0948ba;hpb=9e1051ad0a532338a6f995b9f41f17ac5cdc47a6 diff --git a/ginac/mul.cpp b/ginac/mul.cpp index 9ed27c6c..16bbc5b2 100644 --- a/ginac/mul.cpp +++ b/ginac/mul.cpp @@ -3,7 +3,7 @@ * Implementation of GiNaC's products of expressions. */ /* - * GiNaC Copyright (C) 1999-2007 Johannes Gutenberg University Mainz, Germany + * GiNaC Copyright (C) 1999-2019 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 @@ -20,11 +20,6 @@ * 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" @@ -34,8 +29,14 @@ #include "lst.h" #include "archive.h" #include "utils.h" +#include "symbol.h" #include "compiler.h" +#include +#include +#include +#include + namespace GiNaC { GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(mul, expairseq, @@ -52,7 +53,6 @@ GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(mul, expairseq, mul::mul() { - tinfo_key = &mul::tinfo_static; } ////////// @@ -63,7 +63,6 @@ mul::mul() mul::mul(const ex & lh, const ex & rh) { - tinfo_key = &mul::tinfo_static; overall_coeff = _ex1; construct_from_2_ex(lh,rh); GINAC_ASSERT(is_canonical()); @@ -71,7 +70,6 @@ mul::mul(const ex & lh, const ex & rh) mul::mul(const exvector & v) { - tinfo_key = &mul::tinfo_static; overall_coeff = _ex1; construct_from_exvector(v); GINAC_ASSERT(is_canonical()); @@ -79,7 +77,6 @@ mul::mul(const exvector & v) mul::mul(const epvector & v) { - tinfo_key = &mul::tinfo_static; overall_coeff = _ex1; construct_from_epvector(v); GINAC_ASSERT(is_canonical()); @@ -87,24 +84,27 @@ mul::mul(const epvector & v) mul::mul(const epvector & v, const ex & oc, bool do_index_renaming) { - tinfo_key = &mul::tinfo_static; overall_coeff = oc; construct_from_epvector(v, do_index_renaming); GINAC_ASSERT(is_canonical()); } -mul::mul(std::auto_ptr vp, const ex & oc, bool do_index_renaming) +mul::mul(epvector && vp) +{ + overall_coeff = _ex1; + construct_from_epvector(std::move(vp)); + GINAC_ASSERT(is_canonical()); +} + +mul::mul(epvector && vp, const ex & oc, bool do_index_renaming) { - tinfo_key = &mul::tinfo_static; - GINAC_ASSERT(vp.get()!=0); overall_coeff = oc; - construct_from_epvector(*vp, do_index_renaming); + construct_from_epvector(std::move(vp), do_index_renaming); GINAC_ASSERT(is_canonical()); } mul::mul(const ex & lh, const ex & mh, const ex & rh) { - tinfo_key = &mul::tinfo_static; exvector factors; factors.reserve(3); factors.push_back(lh); @@ -119,8 +119,6 @@ mul::mul(const ex & lh, const ex & mh, const ex & rh) // archiving ////////// -DEFAULT_ARCHIVING(mul) - ////////// // functions overriding virtual functions from base classes ////////// @@ -154,15 +152,13 @@ void mul::do_print(const print_context & c, unsigned level) const print_overall_coeff(c, "*"); - epvector::const_iterator it = seq.begin(), itend = seq.end(); bool first = true; - while (it != itend) { + for (auto & it : seq) { if (!first) c.s << '*'; else first = false; - recombine_pair_to_ex(*it).print(c, precedence()); - ++it; + recombine_pair_to_ex(it).print(c, precedence()); } if (precedence() <= level) @@ -178,15 +174,13 @@ void mul::do_print_latex(const print_latex & c, unsigned level) const // Separate factors into those with negative numeric exponent // and all others - epvector::const_iterator it = seq.begin(), itend = seq.end(); exvector neg_powers, others; - while (it != itend) { - GINAC_ASSERT(is_exactly_a(it->coeff)); - if (ex_to(it->coeff).is_negative()) - neg_powers.push_back(recombine_pair_to_ex(expair(it->rest, -(it->coeff)))); + for (auto & it : seq) { + GINAC_ASSERT(is_exactly_a(it.coeff)); + if (ex_to(it.coeff).is_negative()) + neg_powers.push_back(recombine_pair_to_ex(expair(it.rest, -it.coeff))); else - others.push_back(recombine_pair_to_ex(*it)); - ++it; + others.push_back(recombine_pair_to_ex(it)); } if (!neg_powers.empty()) { @@ -201,11 +195,9 @@ void mul::do_print_latex(const print_latex & c, unsigned level) const } else { // All other factors are printed in the ordinary way - exvector::const_iterator vit = others.begin(), vitend = others.end(); - while (vit != vitend) { + for (auto & vit : others) { c.s << ' '; - vit->print(c, precedence()); - ++vit; + vit.print(c, precedence()); } } @@ -228,7 +220,7 @@ void mul::do_print_csrc(const print_csrc & c, unsigned level) const } // Print arguments, separated by "*" or "/" - epvector::const_iterator it = seq.begin(), itend = seq.end(); + auto it = seq.begin(), itend = seq.end(); while (it != itend) { // If the first argument is a negative integer power, it gets printed as "1.0/" @@ -245,11 +237,9 @@ void mul::do_print_csrc(const print_csrc & c, unsigned level) const if (it->coeff.is_equal(_ex1) || it->coeff.is_equal(_ex_1)) it->rest.print(c, precedence()); else if (it->coeff.info(info_flags::negint)) - // Outer parens around ex needed for broken GCC parser: - (ex(power(it->rest, -ex_to(it->coeff)))).print(c, level); + ex(power(it->rest, -ex_to(it->coeff))).print(c, level); else - // Outer parens around ex needed for broken GCC parser: - (ex(power(it->rest, ex_to(it->coeff)))).print(c, level); + ex(power(it->rest, ex_to(it->coeff))).print(c, level); if (needclosingparenthesis) c.s << ")"; @@ -286,42 +276,135 @@ 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::even: case info_flags::crational_polynomial: case info_flags::rational_function: { - epvector::const_iterator i = seq.begin(), end = seq.end(); - while (i != end) { - if (!(recombine_pair_to_ex(*i).info(inf))) + for (auto & it : seq) { + if (!recombine_pair_to_ex(it).info(inf)) 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: { - epvector::const_iterator i = seq.begin(), end = seq.end(); - while (i != end) { - if ((recombine_pair_to_ex(*i).info(inf))) - return true; - ++i; + case info_flags::positive: + case info_flags::negative: { + if ((inf==info_flags::positive) && (flags & status_flags::is_positive)) + return true; + else if ((inf==info_flags::negative) && (flags & status_flags::is_negative)) + return true; + if (flags & status_flags::purely_indefinite) + return false; + + bool pos = true; + for (auto & it : seq) { + const ex& factor = recombine_pair_to_ex(it); + if (factor.info(info_flags::positive)) + continue; + else if (factor.info(info_flags::negative)) + pos = !pos; + else + return false; } - return false; + if (overall_coeff.info(info_flags::negative)) + pos = !pos; + setflag(pos ? status_flags::is_positive : status_flags::is_negative); + return (inf == info_flags::positive? pos : !pos); + } + case info_flags::nonnegative: { + if (flags & status_flags::is_positive) + return true; + bool pos = true; + for (auto & it : seq) { + const ex& factor = recombine_pair_to_ex(it); + if (factor.info(info_flags::nonnegative) || factor.info(info_flags::positive)) + continue; + else if (factor.info(info_flags::negative)) + pos = !pos; + else + return false; + } + return (overall_coeff.info(info_flags::negative)? !pos : pos); + } + case info_flags::posint: + case info_flags::negint: { + bool pos = true; + for (auto & it : seq) { + const ex& factor = recombine_pair_to_ex(it); + if (factor.info(info_flags::posint)) + continue; + else if (factor.info(info_flags::negint)) + pos = !pos; + else + return false; + } + if (overall_coeff.info(info_flags::negint)) + pos = !pos; + else if (!overall_coeff.info(info_flags::posint)) + return false; + return (inf ==info_flags::posint? pos : !pos); + } + case info_flags::nonnegint: { + bool pos = true; + for (auto & it : seq) { + const ex& factor = recombine_pair_to_ex(it); + if (factor.info(info_flags::nonnegint) || factor.info(info_flags::posint)) + continue; + else if (factor.info(info_flags::negint)) + pos = !pos; + else + return false; + } + if (overall_coeff.info(info_flags::negint)) + pos = !pos; + else if (!overall_coeff.info(info_flags::posint)) + return false; + return pos; + } + case info_flags::indefinite: { + if (flags & status_flags::purely_indefinite) + return true; + if (flags & (status_flags::is_positive | status_flags::is_negative)) + return false; + for (auto & it : seq) { + const ex& term = recombine_pair_to_ex(it); + if (term.info(info_flags::positive) || term.info(info_flags::negative)) + return false; + } + setflag(status_flags::purely_indefinite); + return true; } } return inherited::info(inf); } +bool mul::is_polynomial(const ex & var) const +{ + for (auto & it : seq) { + if (!it.rest.is_polynomial(var) || + (it.rest.has(var) && !it.coeff.info(info_flags::nonnegint))) { + return false; + } + } + return true; +} + int mul::degree(const ex & s) const { // Sum up degrees of factors int deg_sum = 0; - epvector::const_iterator i = seq.begin(), end = seq.end(); - while (i != end) { - if (ex_to(i->coeff).is_integer()) - deg_sum += recombine_pair_to_ex(*i).degree(s); + for (auto & it : seq) { + if (ex_to(it.coeff).is_integer()) + deg_sum += recombine_pair_to_ex(it).degree(s); else { - if (i->rest.has(s)) + if (it.rest.has(s)) throw std::runtime_error("mul::degree() undefined degree because of non-integer exponent"); } - ++i; } return deg_sum; } @@ -330,15 +413,13 @@ int mul::ldegree(const ex & s) const { // Sum up degrees of factors int deg_sum = 0; - epvector::const_iterator i = seq.begin(), end = seq.end(); - while (i != end) { - if (ex_to(i->coeff).is_integer()) - deg_sum += recombine_pair_to_ex(*i).ldegree(s); + for (auto & it : seq) { + if (ex_to(it.coeff).is_integer()) + deg_sum += recombine_pair_to_ex(it).ldegree(s); else { - if (i->rest.has(s)) + if (it.rest.has(s)) throw std::runtime_error("mul::ldegree() undefined degree because of non-integer exponent"); } - ++i; } return deg_sum; } @@ -351,19 +432,15 @@ ex mul::coeff(const ex & s, int n) const if (n==0) { // product of individual coeffs // if a non-zero power of s is found, the resulting product will be 0 - epvector::const_iterator i = seq.begin(), end = seq.end(); - while (i != end) { - coeffseq.push_back(recombine_pair_to_ex(*i).coeff(s,n)); - ++i; - } + for (auto & it : seq) + coeffseq.push_back(recombine_pair_to_ex(it).coeff(s,n)); coeffseq.push_back(overall_coeff); - return (new mul(coeffseq))->setflag(status_flags::dynallocated); + return dynallocate(coeffseq); } - epvector::const_iterator i = seq.begin(), end = seq.end(); bool coeff_found = false; - while (i != end) { - ex t = recombine_pair_to_ex(*i); + for (auto & it : seq) { + ex t = recombine_pair_to_ex(it); ex c = t.coeff(s, n); if (!c.is_zero()) { coeffseq.push_back(c); @@ -371,11 +448,10 @@ ex mul::coeff(const ex & s, int n) const } else { coeffseq.push_back(t); } - ++i; } if (coeff_found) { coeffseq.push_back(overall_coeff); - return (new mul(coeffseq))->setflag(status_flags::dynallocated); + return dynallocate(coeffseq); } return _ex0; @@ -388,41 +464,21 @@ ex mul::coeff(const ex & s, int n) const * - *(+(x1,x2,...);c) -> *(+(*(x1,c),*(x2,c),...)) * - *(x;1) -> x * - *(;c) -> c - * - * @param level cut-off in recursive evaluation */ -ex mul::eval(int level) const + */ +ex mul::eval() const { - std::auto_ptr evaled_seqp = evalchildren(level); - if (evaled_seqp.get()) { - // do more evaluation later - return (new mul(evaled_seqp, overall_coeff))-> - setflag(status_flags::dynallocated); - } - -#ifdef DO_GINAC_ASSERT - epvector::const_iterator i = seq.begin(), end = seq.end(); - while (i != end) { - GINAC_ASSERT((!is_exactly_a(i->rest)) || - (!(ex_to(i->coeff).is_integer()))); - GINAC_ASSERT(!(i->is_canonical_numeric())); - if (is_exactly_a(recombine_pair_to_ex(*i))) - print(print_tree(std::cerr)); - GINAC_ASSERT(!is_exactly_a(recombine_pair_to_ex(*i))); - /* for paranoia */ - expair p = split_ex_to_pair(recombine_pair_to_ex(*i)); - GINAC_ASSERT(p.rest.is_equal(i->rest)); - GINAC_ASSERT(p.coeff.is_equal(i->coeff)); - /* end paranoia */ - ++i; - } -#endif // def DO_GINAC_ASSERT - if (flags & status_flags::evaluated) { GINAC_ASSERT(seq.size()>0); GINAC_ASSERT(seq.size()>1 || !overall_coeff.is_equal(_ex1)); return *this; } - + + const epvector evaled = evalchildren(); + if (unlikely(!evaled.empty())) { + // start over evaluating a new object + return dynallocate(std::move(evaled), overall_coeff); + } + size_t seq_size = seq.size(); if (overall_coeff.is_zero()) { // *(...,x;0) -> 0 @@ -438,25 +494,21 @@ ex mul::eval(int level) const 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); - distrseq->reserve(addref.seq.size()); - epvector::const_iterator i = addref.seq.begin(), end = addref.seq.end(); - while (i != end) { - distrseq->push_back(addref.combine_pair_with_coeff_to_pair(*i, overall_coeff)); - ++i; + epvector distrseq; + distrseq.reserve(addref.seq.size()); + for (auto & it : addref.seq) { + distrseq.push_back(addref.combine_pair_with_coeff_to_pair(it, overall_coeff)); } - return (new add(distrseq, - ex_to(addref.overall_coeff). - mul_dyn(ex_to(overall_coeff))) - )->setflag(status_flags::dynallocated | status_flags::evaluated); + return dynallocate(std::move(distrseq), + ex_to(addref.overall_coeff).mul_dyn(ex_to(overall_coeff))) + .setflag(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 + // 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); + auto i = seq.begin(), last = seq.end(); + auto j = seq.begin(); + epvector s; numeric oc = *_num1_p; bool something_changed = false; while (i!=last) { @@ -468,8 +520,8 @@ ex mul::eval(int level) const // 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_dyn(c); + 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 @@ -483,12 +535,12 @@ ex mul::eval(int level) const } if (! something_changed) { - s->reserve(seq_size); + s.reserve(seq_size); something_changed = true; } while ((j!=i) && (j!=last)) { - s->push_back(*j); + s.push_back(*j); ++j; } @@ -499,62 +551,48 @@ ex mul::eval(int level) const // 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)); + add & primitive = dynallocate(addref); + primitive.clearflag(status_flags::hash_calculated); + primitive.overall_coeff = ex_to(primitive.overall_coeff).div_dyn(c); + for (auto & ai : primitive.seq) + 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); + s.push_back(*j); ++j; } - return (new mul(s, ex_to(overall_coeff).mul_dyn(oc)) - )->setflag(status_flags::dynallocated); + return dynallocate(std::move(s), ex_to(overall_coeff).mul_dyn(oc)); } } return this->hold(); } -ex mul::evalf(int level) const +ex mul::evalf() const { - if (level==1) - return mul(seq,overall_coeff); - - if (level==-max_recursion_level) - throw(std::runtime_error("max recursion level reached")); - - std::auto_ptr s(new epvector); - s->reserve(seq.size()); + epvector s; + s.reserve(seq.size()); - --level; - epvector::const_iterator i = seq.begin(), end = seq.end(); - while (i != end) { - s->push_back(combine_ex_with_coeff_to_pair(i->rest.evalf(level), - i->coeff)); - ++i; - } - return mul(s, overall_coeff.evalf(level)); + for (auto & it : seq) + s.push_back(expair(it.rest.evalf(), it.coeff)); + return dynallocate(std::move(s), overall_coeff.evalf()); } 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); + for (auto & it : seq) { + ex factor = recombine_pair_to_ex(it); ex new_rp = factor.real_part(); ex new_ip = factor.imag_part(); - if(new_ip.is_zero()) { + if (new_ip.is_zero()) { rp *= new_rp; ip *= new_rp; } else { @@ -591,21 +629,19 @@ ex mul::evalm() const // Evaluate children first, look whether there are any matrices at all // (there can be either no matrices or one matrix; if there were more // than one matrix, it would be a non-commutative product) - std::auto_ptr s(new epvector); - s->reserve(seq.size()); + epvector s; + s.reserve(seq.size()); bool have_matrix = false; epvector::iterator the_matrix; - epvector::const_iterator i = seq.begin(), end = seq.end(); - while (i != end) { - const ex &m = recombine_pair_to_ex(*i).evalm(); - s->push_back(split_ex_to_pair(m)); + for (auto & it : seq) { + const ex &m = recombine_pair_to_ex(it).evalm(); + s.push_back(split_ex_to_pair(m)); if (is_a(m)) { have_matrix = true; - the_matrix = s->end() - 1; + the_matrix = s.end() - 1; } - ++i; } if (have_matrix) { @@ -613,12 +649,12 @@ ex mul::evalm() const // The product contained a matrix. We will multiply all other factors // into that matrix. matrix m = ex_to(the_matrix->rest); - s->erase(the_matrix); - ex scalar = (new mul(s, overall_coeff))->setflag(status_flags::dynallocated); + s.erase(the_matrix); + ex scalar = dynallocate(std::move(s), overall_coeff); return m.mul_scalar(scalar); } else - return (new mul(s, overall_coeff))->setflag(status_flags::dynallocated); + return dynallocate(std::move(s), overall_coeff); } ex mul::eval_ncmul(const exvector & v) const @@ -627,16 +663,13 @@ ex mul::eval_ncmul(const exvector & v) const return inherited::eval_ncmul(v); // Find first noncommutative element and call its eval_ncmul() - epvector::const_iterator i = seq.begin(), end = seq.end(); - while (i != end) { - if (i->rest.return_type() == return_types::noncommutative) - return i->rest.eval_ncmul(v); - ++i; - } + for (auto & it : seq) + if (it.rest.return_type() == return_types::noncommutative) + return it.rest.eval_ncmul(v); 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; @@ -668,7 +701,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; @@ -679,7 +712,7 @@ 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) +/** Checks whether 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 @@ -687,17 +720,20 @@ bool tryfactsubs(const ex & origfactor, const ex & patternfactor, int & nummatch * 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, lst &repls, - int factor, int &nummatches, const std::vector &subsed, - std::vector &matched) +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) { - if (factor == pat.nops()) + 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)) { - lst repls; + exmap repls; int nummatches = std::numeric_limits::max(); - std::vector subsed(seq.size(), false); - std::vector matched(seq.size(), false); + 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; @@ -733,46 +769,45 @@ bool mul::has(const ex & pattern, unsigned options) const 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) { + for (auto & it : m) { - if (is_exactly_a(it->first)) { + if (is_exactly_a(it.first)) { retry1: int nummatches = std::numeric_limits::max(); - std::vector currsubsed(seq.size(), false); - lst repls; + std::vector currsubsed(nops(), false); + exmap repls; - if(!algebraic_match_mul_with_mul(*this, it->first, repls, 0, nummatches, subsed, currsubsed)) + if (!algebraic_match_mul_with_mul(*this, it.first, repls, 0, nummatches, subsed, currsubsed)) continue; for (size_t j=0; jfirst.subs(ex(repls), subs_options::no_pattern); - divide_by *= power(subsed_pattern, nummatches); + = it.first.subs(repls, subs_options::no_pattern); + divide_by *= pow(subsed_pattern, nummatches); ex subsed_result - = it->second.subs(ex(repls), subs_options::no_pattern); - multiply_by *= power(subsed_result, nummatches); + = it.second.subs(repls, subs_options::no_pattern); + multiply_by *= pow(subsed_result, nummatches); goto retry1; } else { for (size_t j=0; jnops(); j++) { int nummatches = std::numeric_limits::max(); - lst repls; - if (!subsed[j] && tryfactsubs(op(j), it->first, nummatches, repls)){ + exmap repls; + if (!subsed[j] && tryfactsubs(op(j), it.first, nummatches, repls)){ subsed[j] = true; ex subsed_pattern - = it->first.subs(ex(repls), subs_options::no_pattern); - divide_by *= power(subsed_pattern, nummatches); + = it.first.subs(repls, subs_options::no_pattern); + divide_by *= pow(subsed_pattern, nummatches); ex subsed_result - = it->second.subs(ex(repls), subs_options::no_pattern); - multiply_by *= power(subsed_result, nummatches); + = it.second.subs(repls, subs_options::no_pattern); + multiply_by *= pow(subsed_result, nummatches); } } } @@ -791,6 +826,36 @@ retry1: return ((*this)/divide_by)*multiply_by; } +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. + std::unique_ptr newepv(nullptr); + for (auto 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.reset(new epvector); + newepv->reserve(seq.size()); + for (auto 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; + } + return thisexpairseq(newepv ? std::move(*newepv) : seq, x); +} + + // protected /** Implementation of ex::diff() for a product. It applies the product rule. @@ -803,17 +868,17 @@ ex mul::derivative(const symbol & s) const // D(a*b*c) = D(a)*b*c + a*D(b)*c + a*b*D(c) epvector mulseq = seq; - epvector::const_iterator i = seq.begin(), end = seq.end(); - epvector::iterator i2 = mulseq.begin(); + auto i = seq.begin(), end = seq.end(); + auto i2 = mulseq.begin(); while (i != end) { - expair ep = split_ex_to_pair(power(i->rest, i->coeff - _ex1) * + expair ep = split_ex_to_pair(pow(i->rest, i->coeff - _ex1) * i->rest.diff(s)); ep.swap(*i2); - addseq.push_back((new mul(mulseq, overall_coeff * i->coeff))->setflag(status_flags::dynallocated)); + addseq.push_back(dynallocate(mulseq, overall_coeff * i->coeff)); ep.swap(*i2); ++i; ++i2; } - return (new add(addseq))->setflag(status_flags::dynallocated); + return dynallocate(addseq); } int mul::compare_same_type(const basic & other) const @@ -853,31 +918,29 @@ unsigned mul::return_type() const // all factors checked return all_commutative ? return_types::commutative : return_types::noncommutative; } - -tinfo_t mul::return_type_tinfo() const + +return_type_t mul::return_type_tinfo() const { if (seq.empty()) - return this; // 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(); - while (i != end) { - if (i->rest.return_type() == return_types::noncommutative) - return i->rest.return_type_tinfo(); - ++i; - } + for (auto & it : seq) + if (it.rest.return_type() == return_types::noncommutative) + return it.rest.return_type_tinfo(); + // no noncommutative element found, should not happen - return this; + return make_return_type_t(); } ex mul::thisexpairseq(const epvector & v, const ex & oc, bool do_index_renaming) const { - return (new mul(v, oc, do_index_renaming))->setflag(status_flags::dynallocated); + return dynallocate(v, oc, do_index_renaming); } -ex mul::thisexpairseq(std::auto_ptr vp, const ex & oc, bool do_index_renaming) const +ex mul::thisexpairseq(epvector && vp, const ex & oc, bool do_index_renaming) const { - return (new mul(vp, oc, do_index_renaming))->setflag(status_flags::dynallocated); + return dynallocate(std::move(vp), oc, do_index_renaming); } expair mul::split_ex_to_pair(const ex & e) const @@ -889,60 +952,77 @@ expair mul::split_ex_to_pair(const ex & e) const } return expair(e,_ex1); } - + expair mul::combine_ex_with_coeff_to_pair(const ex & e, const ex & c) const { + GINAC_ASSERT(is_exactly_a(c)); + + // First, try a common shortcut: + if (is_exactly_a(e)) + return expair(e, c); + + // trivial case: exponent 1 + if (c.is_equal(_ex1)) + return split_ex_to_pair(e); + // to avoid duplication of power simplification rules, // we create a temporary power object // otherwise it would be hard to correctly evaluate // expression like (4^(1/3))^(3/2) - if (c.is_equal(_ex1)) - return split_ex_to_pair(e); - - return split_ex_to_pair(power(e,c)); + return split_ex_to_pair(pow(e,c)); } - + expair mul::combine_pair_with_coeff_to_pair(const expair & p, const ex & c) const { + GINAC_ASSERT(is_exactly_a(p.coeff)); + GINAC_ASSERT(is_exactly_a(c)); + + // First, try a common shortcut: + if (is_exactly_a(p.rest)) + return expair(p.rest, p.coeff * c); + + // trivial case: exponent 1 + if (c.is_equal(_ex1)) + return p; + if (p.coeff.is_equal(_ex1)) + return expair(p.rest, c); + // to avoid duplication of power simplification rules, // we create a temporary power object // otherwise it would be hard to correctly evaluate // expression like (4^(1/3))^(3/2) - if (c.is_equal(_ex1)) - return p; - - return split_ex_to_pair(power(recombine_pair_to_ex(p),c)); + return split_ex_to_pair(pow(recombine_pair_to_ex(p),c)); } - + ex mul::recombine_pair_to_ex(const expair & p) const { - if (ex_to(p.coeff).is_equal(*_num1_p)) + if (p.coeff.is_equal(_ex1)) return p.rest; else - return (new power(p.rest,p.coeff))->setflag(status_flags::dynallocated); + return dynallocate(p.rest, p.coeff); } bool mul::expair_needs_further_processing(epp it) { if (is_exactly_a(it->rest) && - ex_to(it->coeff).is_integer()) { + ex_to(it->coeff).is_integer()) { // combined pair is product with integer power -> expand it *it = split_ex_to_pair(recombine_pair_to_ex(*it)); return true; } if (is_exactly_a(it->rest)) { + if (it->coeff.is_equal(_ex1)) { + // pair has coeff 1 and must be moved to the end + return true; + } expair ep = split_ex_to_pair(recombine_pair_to_ex(*it)); if (!ep.is_equal(*it)) { // combined pair is a numeric power which can be simplified *it = ep; return true; } - if (it->coeff.is_equal(_ex1)) { - // combined pair has coeff 1 and must be moved to the end - return true; - } } return false; } @@ -970,17 +1050,16 @@ void mul::combine_overall_coeff(const ex & c1, const ex & c2) bool mul::can_make_flat(const expair & p) const { GINAC_ASSERT(is_exactly_a(p.coeff)); - // 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_p); + + // (x*y)^c == x^c*y^c if c ∈ ℤ + return p.coeff.info(info_flags::integer); } bool mul::can_be_further_expanded(const ex & e) { if (is_exactly_a(e)) { - for (epvector::const_iterator cit = ex_to(e).seq.begin(); cit != ex_to(e).seq.end(); ++cit) { - if (is_exactly_a(cit->rest) && cit->coeff.info(info_flags::posint)) + for (auto & it : ex_to(e).seq) { + if (is_exactly_a(it.rest) && it.coeff.info(info_flags::posint)) return true; } } else if (is_exactly_a(e)) { @@ -992,9 +1071,29 @@ bool mul::can_be_further_expanded(const ex & e) ex mul::expand(unsigned options) const { + // Check for trivial case: expanding the monomial (~ 30% of all calls) + bool monomial_case = true; + for (const auto & i : seq) { + if (!is_a(i.rest) || !i.coeff.info(info_flags::integer)) { + monomial_case = false; + break; + } + } + if (monomial_case) { + 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); + epvector expanded = expandchildren(options); + const epvector & expanded_seq = (expanded.empty() ? seq : expanded); // Now, look for all the factors that are sums and multiply each one out // with the next one that is found while collecting the factors which are @@ -1004,93 +1103,94 @@ ex mul::expand(unsigned options) const epvector non_adds; non_adds.reserve(expanded_seq.size()); - for (epvector::const_iterator cit = expanded_seq.begin(); cit != expanded_seq.end(); ++cit) { - if (is_exactly_a(cit->rest) && - (cit->coeff.is_equal(_ex1))) { + for (const auto & cit : expanded_seq) { + if (is_exactly_a(cit.rest) && + (cit.coeff.is_equal(_ex1))) { if (is_exactly_a(last_expanded)) { // Expand a product of two sums, aggressive version. // Caring for the overall coefficients in separate loops can // sometimes give a performance gain of up to 15%! - const int sizedifference = ex_to(last_expanded).seq.size()-ex_to(cit->rest).seq.size(); + const int sizedifference = ex_to(last_expanded).seq.size()-ex_to(cit.rest).seq.size(); // add2 is for the inner loop and should be the bigger of the two sums // in the presence of asymptotically good sorting: - const add& add1 = (sizedifference<0 ? ex_to(last_expanded) : ex_to(cit->rest)); - const add& add2 = (sizedifference<0 ? ex_to(cit->rest) : ex_to(last_expanded)); - const epvector::const_iterator add1begin = add1.seq.begin(); - const epvector::const_iterator add1end = add1.seq.end(); - const epvector::const_iterator add2begin = add2.seq.begin(); - const epvector::const_iterator add2end = add2.seq.end(); + const add& add1 = (sizedifference<0 ? ex_to(last_expanded) : ex_to(cit.rest)); + const add& add2 = (sizedifference<0 ? ex_to(cit.rest) : ex_to(last_expanded)); epvector distrseq; distrseq.reserve(add1.seq.size()+add2.seq.size()); // Multiply add2 with the overall coefficient of add1 and append it to distrseq: if (!add1.overall_coeff.is_zero()) { if (add1.overall_coeff.is_equal(_ex1)) - distrseq.insert(distrseq.end(),add2begin,add2end); + distrseq.insert(distrseq.end(), add2.seq.begin(), add2.seq.end()); else - for (epvector::const_iterator i=add2begin; i!=add2end; ++i) - distrseq.push_back(expair(i->rest, ex_to(i->coeff).mul_dyn(ex_to(add1.overall_coeff)))); + for (const auto & i : add2.seq) + distrseq.push_back(expair(i.rest, ex_to(i.coeff).mul_dyn(ex_to(add1.overall_coeff)))); } // Multiply add1 with the overall coefficient of add2 and append it to distrseq: if (!add2.overall_coeff.is_zero()) { if (add2.overall_coeff.is_equal(_ex1)) - distrseq.insert(distrseq.end(),add1begin,add1end); + distrseq.insert(distrseq.end(), add1.seq.begin(), add1.seq.end()); else - for (epvector::const_iterator i=add1begin; i!=add1end; ++i) - distrseq.push_back(expair(i->rest, ex_to(i->coeff).mul_dyn(ex_to(add2.overall_coeff)))); + for (const auto & i : add1.seq) + distrseq.push_back(expair(i.rest, ex_to(i.coeff).mul_dyn(ex_to(add2.overall_coeff)))); } // 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); + ex tmp_accu = dynallocate(distrseq, add1.overall_coeff*add2.overall_coeff); 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 (const auto & i : add1.seq) { + add_indices = get_all_dummy_indices_safely(i.rest); + add1_dummy_indices.insert(add1_dummy_indices.end(), add_indices.begin(), add_indices.end()); + } + for (const auto & i : add2.seq) { + 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) { + for (const auto & i2 : add2.seq) { // 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); - for (epvector::const_iterator i1=add1begin; i1!=add1end; ++i1) { + 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 (const auto & i1 : add1.seq) { // 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_new))->setflag(status_flags::dynallocated); + const ex rest = dynallocate(i1.rest, i2_new); if (is_exactly_a(rest)) { - oc += ex_to(rest).mul(ex_to(i1->coeff).mul(ex_to(i2->coeff))); + 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)))); + 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 += dynallocate(std::move(distrseq2), oc); } last_expanded = tmp_accu; - } else { if (!last_expanded.is_equal(_ex1)) non_adds.push_back(split_ex_to_pair(last_expanded)); - last_expanded = cit->rest; + last_expanded = cit.rest; } } else { - non_adds.push_back(*cit); + non_adds.push_back(cit); } } @@ -1100,13 +1200,19 @@ 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; isetflag(status_flags::dynallocated); + 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 = dynallocate(factors, overall_coeff); if (can_be_further_expanded(term)) { distrseq.push_back(term.expand()); } else { @@ -1116,12 +1222,11 @@ ex mul::expand(unsigned options) const } } - return ((new add(distrseq))-> - setflag(status_flags::dynallocated | (options == 0 ? status_flags::expanded : 0))); + return dynallocate(distrseq).setflag(options == 0 ? status_flags::expanded : 0); } non_adds.push_back(split_ex_to_pair(last_expanded)); - ex result = (new mul(non_adds, overall_coeff))->setflag(status_flags::dynallocated); + ex result = dynallocate(non_adds, overall_coeff); if (can_be_further_expanded(result)) { return result.expand(); } else { @@ -1145,46 +1250,47 @@ ex mul::expand(unsigned options) const /** Member-wise expand the expairs representing this sequence. This must be * overridden from expairseq::expandchildren() and done iteratively in order - * to allow for early cancallations and thus safe memory. + * to allow for early cancellations and thus safe memory. * * @see mul::expand() - * @return pointer to epvector containing expanded representation or zero - * pointer, if sequence is unchanged. */ -std::auto_ptr mul::expandchildren(unsigned options) const + * @return epvector containing expanded pairs, empty if no members + * had to be changed. */ +epvector mul::expandchildren(unsigned options) const { - const epvector::const_iterator last = seq.end(); - epvector::const_iterator cit = seq.begin(); + auto cit = seq.begin(), last = seq.end(); while (cit!=last) { const ex & factor = recombine_pair_to_ex(*cit); const ex & expanded_factor = factor.expand(options); if (!are_ex_trivially_equal(factor,expanded_factor)) { // something changed, copy seq, eval and return it - std::auto_ptr s(new epvector); - s->reserve(seq.size()); + epvector s; + s.reserve(seq.size()); // copy parts of seq which are known not to have changed - epvector::const_iterator cit2 = seq.begin(); + auto cit2 = seq.begin(); while (cit2!=cit) { - s->push_back(*cit2); + s.push_back(*cit2); ++cit2; } // copy first changed element - s->push_back(split_ex_to_pair(expanded_factor)); + s.push_back(split_ex_to_pair(expanded_factor)); ++cit2; // copy rest while (cit2!=last) { - s->push_back(split_ex_to_pair(recombine_pair_to_ex(*cit2).expand(options))); + s.push_back(split_ex_to_pair(recombine_pair_to_ex(*cit2).expand(options))); ++cit2; } return s; } ++cit; } - - return std::auto_ptr(0); // nothing has changed + + return epvector(); // nothing has changed } +GINAC_BIND_UNARCHIVER(mul); + } // namespace GiNaC