X-Git-Url: https://ginac.de/ginac.git//ginac.git?a=blobdiff_plain;ds=sidebyside;f=ginac%2Fmul.cpp;h=235478416c42ba64f95e2b34eb7b1eee9a23e7dc;hb=09f37bdbd46f469b3a8a902a43d0f795c41a89bf;hp=6d2e0212bab1c222689fd999caed2bc95be06fda;hpb=c954fb45427efc034382a2a3aa2a783fea2ab0d3;p=ginac.git diff --git a/ginac/mul.cpp b/ginac/mul.cpp index 6d2e0212..23547841 100644 --- a/ginac/mul.cpp +++ b/ginac/mul.cpp @@ -27,6 +27,7 @@ #include "mul.h" #include "add.h" #include "power.h" +#include "operators.h" #include "matrix.h" #include "archive.h" #include "utils.h" @@ -119,7 +120,6 @@ DEFAULT_ARCHIVING(mul) ////////// // public - void mul::print(const print_context & c, unsigned level) const { if (is_a(c)) { @@ -189,7 +189,7 @@ void mul::print(const print_context & c, unsigned level) const bool first = true; // First print the overall numeric coefficient - numeric coeff = ex_to(overall_coeff); + const numeric &coeff = ex_to(overall_coeff); if (coeff.csgn() == -1) c.s << '-'; if (!coeff.is_equal(_num1) && @@ -507,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()) { @@ -586,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)); } @@ -599,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)); } @@ -612,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_exactly_a((*it).rest) && - 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_exactly_a((*it).rest)) { - 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; } @@ -686,20 +681,78 @@ ex mul::expand(unsigned options) const (cit->coeff.is_equal(_ex1))) { ++number_of_adds; 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;