X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=blobdiff_plain;f=ginac%2Fmul.cpp;h=28233b302ebc1b1485b390c055df10586136c6c7;hp=51668803a88371322901cf9c4a35c2f6011d4295;hb=f29f08a2f2f243538d1bedd558991daf49531eb0;hpb=896ad35db238cf5665ca8d519799fa759045b298 diff --git a/ginac/mul.cpp b/ginac/mul.cpp index 51668803..28233b30 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 @@ -141,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 @@ -149,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: @@ -159,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 << "*"; @@ -686,19 +686,78 @@ ex mul::expand(unsigned options) const (cit->coeff.is_equal(_ex1))) { ++number_of_adds; if (is_ex_exactly_of_type(last_expanded, add)) { +#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_ex_exactly_of_type(rest, numeric)) + 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;