X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=blobdiff_plain;f=ginac%2Fmul.cpp;h=7ca04ef36ecb09696bc7b8affe89c60dfa1942c8;hp=cac1ffb7a14e99b2c9a5825d3aac7ac17dd5eb5d;hb=35287d4fef8dc61a10966091ff662eeb9444f87a;hpb=c378e979c7a9b9807941b27adbac690699abe70c diff --git a/ginac/mul.cpp b/ginac/mul.cpp index cac1ffb7..7ca04ef3 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,6 +20,7 @@ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */ +#include #include #include @@ -28,7 +29,6 @@ #include "power.h" #include "matrix.h" #include "archive.h" -#include "debugmsg.h" #include "utils.h" namespace GiNaC { @@ -36,12 +36,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,34 +55,30 @@ 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(); + overall_coeff = _ex1; construct_from_2_ex(lh,rh); GINAC_ASSERT(is_canonical()); } mul::mul(const exvector & v) { - debugmsg("mul ctor from exvector",LOGLEVEL_CONSTRUCT); tinfo_key = TINFO_mul; - overall_coeff = _ex1(); + overall_coeff = _ex1; construct_from_exvector(v); GINAC_ASSERT(is_canonical()); } mul::mul(const epvector & v) { - debugmsg("mul ctor from epvector",LOGLEVEL_CONSTRUCT); tinfo_key = TINFO_mul; - overall_coeff = _ex1(); + overall_coeff = _ex1; construct_from_epvector(v); GINAC_ASSERT(is_canonical()); } 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 +87,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,14 +97,13 @@ 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); factors.push_back(lh); factors.push_back(mh); factors.push_back(rh); - overall_coeff = _ex1(); + overall_coeff = _ex1; construct_from_exvector(factors); GINAC_ASSERT(is_canonical()); } @@ -129,8 +122,6 @@ DEFAULT_ARCHIVING(mul) void mul::print(const print_context & c, unsigned level) const { - debugmsg("mul print", LOGLEVEL_PRINT); - if (is_a(c)) { inherited::print(c, level); @@ -140,8 +131,8 @@ void mul::print(const print_context & c, unsigned level) const if (precedence() <= level) c.s << "("; - if (!overall_coeff.is_equal(_ex1())) { - overall_coeff.bp->print(c, precedence()); + if (!overall_coeff.is_equal(_ex1)) { + overall_coeff.print(c, precedence()); c.s << "*"; } @@ -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); if (coeff.csgn() == -1) c.s << '-'; - if (!coeff.is_equal(_num1()) && - !coeff.is_equal(_num_1())) { + if (!coeff.is_equal(_num1) && + !coeff.is_equal(_num_1)) { if (coeff.is_rational()) { if (coeff.is_negative()) (-coeff).print(c); @@ -327,7 +326,7 @@ ex mul::coeff(const ex & s, int n) const return (new mul(coeffseq))->setflag(status_flags::dynallocated); } - return _ex0(); + return _ex0; } /** Perform automatic term rewriting rules in this class. In the following @@ -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 @@ -353,12 +350,12 @@ ex mul::eval(int level) const #ifdef DO_GINAC_ASSERT epvector::const_iterator i = seq.begin(), end = seq.end(); while (i != end) { - GINAC_ASSERT((!is_ex_exactly_of_type(i->rest, mul)) || + 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)) print(print_tree(std::cerr)); - GINAC_ASSERT(!is_ex_exactly_of_type(recombine_pair_to_ex(*i), numeric)); + 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)); @@ -370,23 +367,23 @@ ex mul::eval(int level) const if (flags & status_flags::evaluated) { GINAC_ASSERT(seq.size()>0); - GINAC_ASSERT(seq.size()>1 || !overall_coeff.is_equal(_ex1())); + GINAC_ASSERT(seq.size()>1 || !overall_coeff.is_equal(_ex1)); return *this; } int seq_size = seq.size(); if (overall_coeff.is_zero()) { // *(...,x;0) -> 0 - return _ex0(); + return _ex0; } else if (seq_size==0) { // *(;c) -> c return overall_coeff; - } else if (seq_size==1 && overall_coeff.is_equal(_ex1())) { + } else if (seq_size==1 && overall_coeff.is_equal(_ex1)) { // *(x;1) -> x return recombine_pair_to_ex(*(seq.begin())); } else if ((seq_size==1) && is_ex_exactly_of_type((*seq.begin()).rest,add) && - ex_to((*seq.begin()).coeff).is_equal(_num1())) { + 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); epvector *distrseq = new epvector(); @@ -428,7 +425,7 @@ ex mul::evalf(int level) const ex mul::evalm(void) const { // numeric*matrix - if (seq.size() == 1 && seq[0].coeff.is_equal(_ex1()) + if (seq.size() == 1 && seq[0].coeff.is_equal(_ex1) && is_ex_of_type(seq[0].rest, matrix)) return ex_to(seq[0].rest).mul(ex_to(overall_coeff)); @@ -495,7 +492,7 @@ ex mul::derivative(const symbol & s) const epvector::const_iterator i = seq.begin(), end = seq.end(); epvector::iterator i2 = mulseq.begin(); while (i != end) { - expair ep = split_ex_to_pair(power(i->rest, i->coeff - _ex1()) * + expair ep = split_ex_to_pair(power(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)); @@ -581,7 +578,7 @@ expair mul::split_ex_to_pair(const ex & e) const if (is_ex_exactly_of_type(powerref.exponent,numeric)) return expair(powerref.basis,powerref.exponent); } - return expair(e,_ex1()); + return expair(e,_ex1); } expair mul::combine_ex_with_coeff_to_pair(const ex & e, @@ -591,7 +588,7 @@ expair mul::combine_ex_with_coeff_to_pair(const ex & e, // we create a temporary power object // otherwise it would be hard to correctly simplify // expression like (4^(1/3))^(3/2) - if (are_ex_trivially_equal(c,_ex1())) + if (are_ex_trivially_equal(c,_ex1)) return split_ex_to_pair(e); return split_ex_to_pair(power(e,c)); @@ -604,7 +601,7 @@ expair mul::combine_pair_with_coeff_to_pair(const expair & p, // we create a temporary power object // otherwise it would be hard to correctly simplify // expression like (4^(1/3))^(3/2) - if (are_ex_trivially_equal(c,_ex1())) + if (are_ex_trivially_equal(c,_ex1)) return p; return split_ex_to_pair(power(recombine_pair_to_ex(p),c)); @@ -612,7 +609,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)) return p.rest; else return power(p.rest,p.coeff); @@ -633,7 +630,7 @@ bool mul::expair_needs_further_processing(epp it) *it = ep; return true; } - if (ex_to((*it).coeff).is_equal(_num1())) { + if (ex_to((*it).coeff).is_equal(_num1)) { // combined pair has coeff 1 and must be moved to the end return true; } @@ -643,31 +640,31 @@ bool mul::expair_needs_further_processing(epp it) ex mul::default_overall_coeff(void) const { - return _ex1(); + return _ex1; } void mul::combine_overall_coeff(const ex & c) { - GINAC_ASSERT(is_ex_exactly_of_type(overall_coeff,numeric)); - GINAC_ASSERT(is_ex_exactly_of_type(c,numeric)); + GINAC_ASSERT(is_exactly_a(overall_coeff)); + GINAC_ASSERT(is_exactly_a(c)); overall_coeff = ex_to(overall_coeff).mul_dyn(ex_to(c)); } void mul::combine_overall_coeff(const ex & c1, const ex & c2) { - GINAC_ASSERT(is_ex_exactly_of_type(overall_coeff,numeric)); - GINAC_ASSERT(is_ex_exactly_of_type(c1,numeric)); - GINAC_ASSERT(is_ex_exactly_of_type(c2,numeric)); + GINAC_ASSERT(is_exactly_a(overall_coeff)); + GINAC_ASSERT(is_exactly_a(c1)); + GINAC_ASSERT(is_exactly_a(c2)); overall_coeff = ex_to(overall_coeff).mul_dyn(ex_to(c1).power(ex_to(c2))); } bool mul::can_make_flat(const expair & p) const { - GINAC_ASSERT(is_ex_exactly_of_type(p.coeff,numeric)); + 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()); + return ex_to(p.coeff).is_equal(_num1); } ex mul::expand(unsigned options) const @@ -680,28 +677,84 @@ ex mul::expand(unsigned options) const // with the next one that is found while collecting the factors which are // not sums int number_of_adds = 0; - ex last_expanded = _ex1(); + ex last_expanded = _ex1; epvector non_adds; 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) && - (cit->coeff.is_equal(_ex1()))) { + (cit->coeff.is_equal(_ex1))) { ++number_of_adds; if (is_ex_exactly_of_type(last_expanded, add)) { +#if 1 + // 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); + distrseq.clear(); + } + last_expanded = tmp_accu; +#else + // Expand a product of two sums, aggressive version. + // Caring for the overall coefficients in separate loops can give + // a performance gain of up to 20%! + const add & add1 = ex_to(last_expanded); + const add & add2 = ex_to(cit->rest); + 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, ex_to(add1.overall_coeff).mul(ex_to(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_ex_exactly_of_type(rest, numeric)) + oc *= ex_to(rest).mul_dyn(ex_to(i1->coeff).mul_dyn(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; @@ -756,7 +809,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);