- exvector sub_expanded_seq;
- intvector positions_of_adds;
- intvector number_of_add_operands;
-
- epvector * expanded_seqp=expandchildren(options);
-
- epvector const & expanded_seq = expanded_seqp==0 ? seq : *expanded_seqp;
-
- positions_of_adds.resize(expanded_seq.size());
- number_of_add_operands.resize(expanded_seq.size());
-
- int number_of_adds=0;
- int number_of_expanded_terms=1;
-
- unsigned current_position=0;
- epvector::const_iterator last=expanded_seq.end();
- for (epvector::const_iterator cit=expanded_seq.begin(); cit!=last; ++cit) {
- if (is_ex_exactly_of_type((*cit).rest,add)&&
- (ex_to_numeric((*cit).coeff).is_equal(_num1()))) {
- positions_of_adds[number_of_adds]=current_position;
- add const & expanded_addref=ex_to_add((*cit).rest);
- unsigned addref_nops=expanded_addref.nops();
- number_of_add_operands[number_of_adds]=addref_nops;
- number_of_expanded_terms *= addref_nops;
- number_of_adds++;
- }
- current_position++;
- }
-
- if (number_of_adds==0) {
- if (expanded_seqp==0) {
- return this->setflag(status_flags::expanded);
- }
- return (new mul(expanded_seqp,overall_coeff))->
- setflag(status_flags::dynallocated ||
- status_flags::expanded);
- }
-
- exvector distrseq;
- distrseq.reserve(number_of_expanded_terms);
-
- intvector k;
- k.resize(number_of_adds);
-
- int l;
- for (l=0; l<number_of_adds; l++) {
- k[l]=0;
- }
-
- while (1) {
- epvector term;
- term=expanded_seq;
- for (l=0; l<number_of_adds; l++) {
- add const & addref=ex_to_add(expanded_seq[positions_of_adds[l]].rest);
- GINAC_ASSERT(term[positions_of_adds[l]].coeff.compare(_ex1())==0);
- term[positions_of_adds[l]]=split_ex_to_pair(addref.op(k[l]));
- }
- /*
- cout << "mul::expand() term begin" << endl;
- for (epvector::const_iterator cit=term.begin(); cit!=term.end(); ++cit) {
- cout << "rest" << endl;
- (*cit).rest.printtree(cout);
- cout << "coeff" << endl;
- (*cit).coeff.printtree(cout);
- }
- cout << "mul::expand() term end" << endl;
- */
- distrseq.push_back((new mul(term,overall_coeff))->
- setflag(status_flags::dynallocated |
- status_flags::expanded));
-
- // increment k[]
- l=number_of_adds-1;
- while ((l>=0)&&((++k[l])>=number_of_add_operands[l])) {
- k[l]=0;
- l--;
- }
- if (l<0) break;
- }
-
- if (expanded_seqp!=0) {
- delete expanded_seqp;
- }
- /*
- cout << "mul::expand() distrseq begin" << endl;
- for (exvector::const_iterator cit=distrseq.begin(); cit!=distrseq.end(); ++cit) {
- (*cit).printtree(cout);
- }
- cout << "mul::expand() distrseq end" << endl;
- */
-
- return (new add(distrseq))->setflag(status_flags::dynallocated |
- status_flags::expanded);
-}
-
+ // First, expand the children
+ epvector * expanded_seqp = expandchildren(options);
+ const epvector & expanded_seq = (expanded_seqp == NULL) ? seq : *expanded_seqp;
+
+ // Now, look for all the factors that are sums and multiply each one out
+ // with the next one that is found while collecting the factors which are
+ // not sums
+ int number_of_adds = 0;
+ 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))) {
+ ++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<add>(last_expanded);
+ const add & add2 = ex_to<add>(cit->rest);
+ const int n1 = add1.nops();
+ const int n2 = add2.nops();
+ ex tmp_accu;
+ exvector distrseq;
+ distrseq.reserve(n2);
+ for (int i1=0; i1<n1; ++i1) {
+ distrseq.clear();
+ // cache the first operand (for efficiency):
+ const ex op1 = add1.op(i1);
+ for (int i2=0; i2<n2; ++i2)
+ distrseq.push_back(op1 * add2.op(i2));
+ tmp_accu += (new add(distrseq))->
+ 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<add>(last_expanded).seq.size()-ex_to<add>(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<add>(last_expanded) : ex_to<add>(cit->rest));
+ const add& add2 = (sizedifference<0 ? ex_to<add>(cit->rest) : ex_to<add>(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<numeric>(i->coeff).mul_dyn(ex_to<numeric>(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<numeric>(i->coeff).mul_dyn(ex_to<numeric>(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<numeric>(rest).mul(ex_to<numeric>(i1->coeff).mul(ex_to<numeric>(i2->coeff)));
+ else
+ distrseq.push_back(expair(rest, ex_to<numeric>(i1->coeff).mul_dyn(ex_to<numeric>(i2->coeff))));
+ }
+ tmp_accu += (new add(distrseq, oc))->setflag(status_flags::dynallocated);
+ }
+ last_expanded = tmp_accu;
+#endif
+ } else {
+ non_adds.push_back(split_ex_to_pair(last_expanded));
+ last_expanded = cit->rest;
+ }
+ } else {
+ non_adds.push_back(*cit);
+ }
+ ++cit;
+ }
+ if (expanded_seqp)
+ delete expanded_seqp;
+
+ // 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)) {
+ const add & finaladd = ex_to<add>(last_expanded);
+ exvector distrseq;
+ int n = finaladd.nops();
+ distrseq.reserve(n);
+ for (int i=0; i<n; ++i) {
+ epvector factors = non_adds;
+ factors.push_back(split_ex_to_pair(finaladd.op(i)));
+ distrseq.push_back((new mul(factors, overall_coeff))->
+ setflag(status_flags::dynallocated | (options == 0 ? status_flags::expanded : 0)));
+ }
+ return ((new add(distrseq))->
+ setflag(status_flags::dynallocated | (options == 0 ? status_flags::expanded : 0)));
+ }
+ non_adds.push_back(split_ex_to_pair(last_expanded));
+ return (new mul(non_adds, overall_coeff))->
+ setflag(status_flags::dynallocated | (options == 0 ? status_flags::expanded : 0));
+}
+
+