X-Git-Url: https://ginac.de/ginac.git//ginac.git?a=blobdiff_plain;f=ginac%2Fmul.cpp;h=235478416c42ba64f95e2b34eb7b1eee9a23e7dc;hb=09f37bdbd46f469b3a8a902a43d0f795c41a89bf;hp=d41179444dce765642e3a8e726586a0e6eb4d992;hpb=27d6204effdef95a00af461fff98024e290dbaa7;p=ginac.git diff --git a/ginac/mul.cpp b/ginac/mul.cpp index d4117944..23547841 100644 --- a/ginac/mul.cpp +++ b/ginac/mul.cpp @@ -3,7 +3,7 @@ * Implementation of GiNaC's products of expressions. */ /* - * GiNaC Copyright (C) 1999-2001 Johannes Gutenberg University Mainz, Germany + * GiNaC Copyright (C) 1999-2002 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,15 +20,16 @@ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */ +#include #include #include #include "mul.h" #include "add.h" #include "power.h" +#include "operators.h" #include "matrix.h" #include "archive.h" -#include "debugmsg.h" #include "utils.h" namespace GiNaC { @@ -36,12 +37,11 @@ namespace GiNaC { GINAC_IMPLEMENT_REGISTERED_CLASS(mul, expairseq) ////////// -// default ctor, dctor, copy ctor assignment operator and helpers +// default ctor, dtor, copy ctor, assignment operator and helpers ////////// mul::mul() { - debugmsg("mul default ctor",LOGLEVEL_CONSTRUCT); tinfo_key = TINFO_mul; } @@ -56,7 +56,6 @@ DEFAULT_DESTROY(mul) mul::mul(const ex & lh, const ex & rh) { - debugmsg("mul ctor from ex,ex",LOGLEVEL_CONSTRUCT); tinfo_key = TINFO_mul; overall_coeff = _ex1; construct_from_2_ex(lh,rh); @@ -65,7 +64,6 @@ mul::mul(const ex & lh, const ex & rh) mul::mul(const exvector & v) { - debugmsg("mul ctor from exvector",LOGLEVEL_CONSTRUCT); tinfo_key = TINFO_mul; overall_coeff = _ex1; construct_from_exvector(v); @@ -74,7 +72,6 @@ mul::mul(const exvector & v) mul::mul(const epvector & v) { - debugmsg("mul ctor from epvector",LOGLEVEL_CONSTRUCT); tinfo_key = TINFO_mul; overall_coeff = _ex1; construct_from_epvector(v); @@ -83,7 +80,6 @@ mul::mul(const epvector & v) mul::mul(const epvector & v, const ex & oc) { - debugmsg("mul ctor from epvector,ex",LOGLEVEL_CONSTRUCT); tinfo_key = TINFO_mul; overall_coeff = oc; construct_from_epvector(v); @@ -92,7 +88,6 @@ mul::mul(const epvector & v, const ex & oc) mul::mul(epvector * vp, const ex & oc) { - debugmsg("mul ctor from epvector *,ex",LOGLEVEL_CONSTRUCT); tinfo_key = TINFO_mul; GINAC_ASSERT(vp!=0); overall_coeff = oc; @@ -103,7 +98,6 @@ mul::mul(epvector * vp, const ex & oc) mul::mul(const ex & lh, const ex & mh, const ex & rh) { - debugmsg("mul ctor from ex,ex,ex",LOGLEVEL_CONSTRUCT); tinfo_key = TINFO_mul; exvector factors; factors.reserve(3); @@ -126,11 +120,8 @@ DEFAULT_ARCHIVING(mul) ////////// // public - void mul::print(const print_context & c, unsigned level) const { - debugmsg("mul print", LOGLEVEL_PRINT); - if (is_a(c)) { inherited::print(c, level); @@ -150,7 +141,7 @@ void mul::print(const print_context & c, unsigned level) const while (it != itend) { // If the first argument is a negative integer power, it gets printed as "1.0/" - if (it == seq.begin() && ex_to(it->coeff).is_integer() && it->coeff.compare(_num0) < 0) { + if (it == seq.begin() && ex_to(it->coeff).is_integer() && it->coeff.info(info_flags::negative)) { if (is_a(c)) c.s << "recip("; else @@ -158,7 +149,7 @@ void mul::print(const print_context & c, unsigned level) const } // If the exponent is 1 or -1, it is left out - if (it->coeff.compare(_ex1) == 0 || it->coeff.compare(_num_1) == 0) + if (it->coeff.is_equal(_ex1) || it->coeff.is_equal(_ex_1)) it->rest.print(c, precedence()); else { // Outer parens around ex needed for broken gcc-2.95 parser: @@ -168,7 +159,7 @@ void mul::print(const print_context & c, unsigned level) const // Separator is "/" for negative integer powers, "*" otherwise ++it; if (it != itend) { - if (ex_to(it->coeff).is_integer() && it->coeff.compare(_num0) < 0) + if (ex_to(it->coeff).is_integer() && it->coeff.info(info_flags::negative)) c.s << "/"; else c.s << "*"; @@ -178,6 +169,14 @@ void mul::print(const print_context & c, unsigned level) const if (precedence() <= level) c.s << ")"; + } else if (is_a(c)) { + c.s << class_name() << '('; + op(0).print(c); + for (unsigned i=1; i(overall_coeff); + const numeric &coeff = ex_to(overall_coeff); if (coeff.csgn() == -1) c.s << '-'; if (!coeff.is_equal(_num1) && @@ -341,8 +340,6 @@ ex mul::coeff(const ex & s, int n) const * @param level cut-off in recursive evaluation */ ex mul::eval(int level) const { - debugmsg("mul eval",LOGLEVEL_MEMBER_FUNCTION); - epvector *evaled_seqp = evalchildren(level); if (evaled_seqp) { // do more evaluation later @@ -356,7 +353,7 @@ ex mul::eval(int level) const GINAC_ASSERT((!is_exactly_a(i->rest)) || (!(ex_to(i->coeff).is_integer()))); GINAC_ASSERT(!(i->is_canonical_numeric())); - if (is_ex_exactly_of_type(recombine_pair_to_ex(*i), 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 */ @@ -385,7 +382,7 @@ ex mul::eval(int level) const // *(x;1) -> x return recombine_pair_to_ex(*(seq.begin())); } else if ((seq_size==1) && - is_ex_exactly_of_type((*seq.begin()).rest,add) && + is_exactly_a((*seq.begin()).rest) && ex_to((*seq.begin()).coeff).is_equal(_num1)) { // *(+(x,y,...);c) -> +(*(x,c),*(y,c),...) (c numeric(), no powers of +()) const add & addref = ex_to((*seq.begin()).rest); @@ -429,7 +426,7 @@ ex mul::evalm(void) const { // numeric*matrix if (seq.size() == 1 && seq[0].coeff.is_equal(_ex1) - && is_ex_of_type(seq[0].rest, matrix)) + && is_a(seq[0].rest)) return ex_to(seq[0].rest).mul(ex_to(overall_coeff)); // Evaluate children first, look whether there are any matrices at all @@ -445,7 +442,7 @@ ex mul::evalm(void) const while (i != end) { const ex &m = recombine_pair_to_ex(*i).evalm(); s->push_back(split_ex_to_pair(m)); - if (is_ex_of_type(m, matrix)) { + if (is_a(m)) { have_matrix = true; the_matrix = s->end() - 1; } @@ -510,11 +507,6 @@ int mul::compare_same_type(const basic & other) const return inherited::compare_same_type(other); } -bool mul::is_equal_same_type(const basic & other) const -{ - return inherited::is_equal_same_type(other); -} - unsigned mul::return_type(void) const { if (seq.empty()) { @@ -576,9 +568,9 @@ ex mul::thisexpairseq(epvector * vp, const ex & oc) const expair mul::split_ex_to_pair(const ex & e) const { - if (is_ex_exactly_of_type(e,power)) { + if (is_exactly_a(e)) { const power & powerref = ex_to(e); - if (is_ex_exactly_of_type(powerref.exponent,numeric)) + if (is_exactly_a(powerref.exponent)) return expair(powerref.basis,powerref.exponent); } return expair(e,_ex1); @@ -589,11 +581,11 @@ expair mul::combine_ex_with_coeff_to_pair(const ex & e, { // to avoid duplication of power simplification rules, // we create a temporary power object - // otherwise it would be hard to correctly simplify + // otherwise it would be hard to correctly evaluate // expression like (4^(1/3))^(3/2) - if (are_ex_trivially_equal(c,_ex1)) + if (c.is_equal(_ex1)) return split_ex_to_pair(e); - + return split_ex_to_pair(power(e,c)); } @@ -602,11 +594,11 @@ expair mul::combine_pair_with_coeff_to_pair(const expair & p, { // to avoid duplication of power simplification rules, // we create a temporary power object - // otherwise it would be hard to correctly simplify + // otherwise it would be hard to correctly evaluate // expression like (4^(1/3))^(3/2) - if (are_ex_trivially_equal(c,_ex1)) + if (c.is_equal(_ex1)) return p; - + return split_ex_to_pair(power(recombine_pair_to_ex(p),c)); } @@ -615,25 +607,25 @@ ex mul::recombine_pair_to_ex(const expair & p) const if (ex_to(p.coeff).is_equal(_num1)) return p.rest; else - return power(p.rest,p.coeff); + return (new power(p.rest,p.coeff))->setflag(status_flags::dynallocated); } bool mul::expair_needs_further_processing(epp it) { - if (is_ex_exactly_of_type((*it).rest,mul) && - ex_to((*it).coeff).is_integer()) { + if (is_exactly_a(it->rest) && + 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_ex_exactly_of_type((*it).rest,numeric)) { - expair ep=split_ex_to_pair(recombine_pair_to_ex(*it)); + if (is_exactly_a(it->rest)) { + 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 (ex_to((*it).coeff).is_equal(_num1)) { + if (it->coeff.is_equal(_ex1)) { // combined pair has coeff 1 and must be moved to the end return true; } @@ -685,23 +677,82 @@ ex mul::expand(unsigned options) const non_adds.reserve(expanded_seq.size()); epvector::const_iterator cit = expanded_seq.begin(), last = expanded_seq.end(); while (cit != last) { - if (is_ex_exactly_of_type(cit->rest, add) && + if (is_exactly_a(cit->rest) && (cit->coeff.is_equal(_ex1))) { ++number_of_adds; - if (is_ex_exactly_of_type(last_expanded, add)) { + if (is_exactly_a(last_expanded)) { +#if 0 + // Expand a product of two sums, simple and robust version. const add & add1 = ex_to(last_expanded); const add & add2 = ex_to(cit->rest); - int n1 = add1.nops(); - int n2 = add2.nops(); + const int n1 = add1.nops(); + const int n2 = add2.nops(); + ex tmp_accu; exvector distrseq; - distrseq.reserve(n1*n2); + distrseq.reserve(n2); for (int i1=0; i1 + setflag(status_flags::dynallocated); + } + last_expanded = tmp_accu; +#else + // 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(); + // 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(); + 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); + 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)))); + } + // 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); + 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)))); + } + // 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); + // Multiply explicitly all non-numeric terms of add1 and add2: + for (epvector::const_iterator i1=add1begin; i1!=add1end; ++i1) { + // 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) { + // 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)) + 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)))); } + tmp_accu += (new add(distrseq, oc))->setflag(status_flags::dynallocated); } - last_expanded = (new add(distrseq))-> - setflag(status_flags::dynallocated | (options == 0 ? status_flags::expanded : 0)); + last_expanded = tmp_accu; +#endif } else { non_adds.push_back(split_ex_to_pair(last_expanded)); last_expanded = cit->rest; @@ -716,7 +767,7 @@ ex mul::expand(unsigned options) const // Now the only remaining thing to do is to multiply the factors which // were not sums into the "last_expanded" sum - if (is_ex_exactly_of_type(last_expanded, add)) { + if (is_exactly_a(last_expanded)) { const add & finaladd = ex_to(last_expanded); exvector distrseq; int n = finaladd.nops(); @@ -756,7 +807,7 @@ ex mul::expand(unsigned options) const * pointer, if sequence is unchanged. */ epvector * mul::expandchildren(unsigned options) const { - epvector::const_iterator last = seq.end(); + const epvector::const_iterator last = seq.end(); epvector::const_iterator cit = seq.begin(); while (cit!=last) { const ex & factor = recombine_pair_to_ex(*cit);