X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=blobdiff_plain;f=ginac%2Fmul.cpp;h=cee5cd7a7a9f4faf3405eb1f4c9af0b52bd26c02;hp=9ed27c6c4bf3f4c7435674de9138bb4aa16d76fc;hb=8ed95601ab0cff1bb02b2a908e5a2c118b1f0a06;hpb=9e1051ad0a532338a6f995b9f41f17ac5cdc47a6;ds=sidebyside diff --git a/ginac/mul.cpp b/ginac/mul.cpp index 9ed27c6c..cee5cd7a 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-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 @@ -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,7 +84,6 @@ 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()); @@ -95,7 +91,6 @@ mul::mul(const epvector & v, const ex & oc, bool do_index_renaming) mul::mul(std::auto_ptr 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); @@ -104,7 +99,6 @@ mul::mul(std::auto_ptr vp, const ex & oc, bool do_index_renaming) 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 +113,6 @@ mul::mul(const ex & lh, const ex & mh, const ex & rh) // archiving ////////// -DEFAULT_ARCHIVING(mul) - ////////// // functions overriding virtual functions from base classes ////////// @@ -286,6 +278,12 @@ 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(); @@ -294,6 +292,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: { @@ -305,10 +305,89 @@ bool mul::info(unsigned inf) const } return false; } + case info_flags::positive: + case info_flags::negative: { + bool pos = true; + 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)) + pos = !pos; + else + return false; + } + if (overall_coeff.info(info_flags::negative)) + pos = !pos; + return (inf ==info_flags::positive? pos : !pos); + } + case info_flags::nonnegative: { + bool pos = true; + 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::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; + 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)) + 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; + 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::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; + } } 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 @@ -399,24 +478,6 @@ ex mul::eval(int level) const 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)); @@ -451,7 +512,7 @@ ex mul::eval(int level) const )->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 + // 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(); @@ -468,8 +529,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 @@ -503,8 +564,7 @@ ex mul::eval(int level) const 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) + 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)); @@ -636,7 +696,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; @@ -668,7 +728,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; @@ -687,17 +747,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,8 +796,7 @@ 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; @@ -743,8 +805,8 @@ ex mul::algebraic_subs_mul(const exmap & m, unsigned options) const 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)) continue; @@ -753,10 +815,10 @@ retry1: if (currsubsed[j]) subsed[j] = true; ex subsed_pattern - = it->first.subs(ex(repls), subs_options::no_pattern); + = it->first.subs(repls, subs_options::no_pattern); divide_by *= power(subsed_pattern, nummatches); ex subsed_result - = it->second.subs(ex(repls), subs_options::no_pattern); + = it->second.subs(repls, subs_options::no_pattern); multiply_by *= power(subsed_result, nummatches); goto retry1; @@ -764,14 +826,14 @@ retry1: for (size_t j=0; jnops(); j++) { int nummatches = std::numeric_limits::max(); - lst 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); + = it->first.subs(repls, subs_options::no_pattern); divide_by *= power(subsed_pattern, nummatches); ex subsed_result - = it->second.subs(ex(repls), subs_options::no_pattern); + = it->second.subs(repls, subs_options::no_pattern); multiply_by *= power(subsed_result, nummatches); } } @@ -791,6 +853,38 @@ 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. + 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. @@ -853,11 +947,11 @@ 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(); @@ -867,7 +961,7 @@ tinfo_t mul::return_type_tinfo() const ++i; } // 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 @@ -889,7 +983,7 @@ 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 { @@ -902,7 +996,7 @@ expair mul::combine_ex_with_coeff_to_pair(const ex & e, return split_ex_to_pair(power(e,c)); } - + expair mul::combine_pair_with_coeff_to_pair(const expair & p, const ex & c) const { @@ -915,7 +1009,7 @@ expair mul::combine_pair_with_coeff_to_pair(const expair & p, return split_ex_to_pair(power(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)) @@ -927,22 +1021,22 @@ ex mul::recombine_pair_to_ex(const expair & p) const 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; } @@ -992,6 +1086,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); @@ -1047,28 +1159,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(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: @@ -1076,13 +1194,12 @@ ex mul::expand(unsigned options) const 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)))); + 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)); @@ -1100,12 +1217,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; isetflag(status_flags::dynallocated); if (can_be_further_expanded(term)) { distrseq.push_back(term.expand()); @@ -1187,4 +1310,6 @@ std::auto_ptr mul::expandchildren(unsigned options) const return std::auto_ptr(0); // nothing has changed } +GINAC_BIND_UNARCHIVER(mul); + } // namespace GiNaC