* As agreed upon at lunch, remove debugmsg since it is hardly of any use.
[ginac.git] / ginac / mul.cpp
1 /** @file mul.cpp
2  *
3  *  Implementation of GiNaC's products of expressions. */
4
5 /*
6  *  GiNaC Copyright (C) 1999-2001 Johannes Gutenberg University Mainz, Germany
7  *
8  *  This program is free software; you can redistribute it and/or modify
9  *  it under the terms of the GNU General Public License as published by
10  *  the Free Software Foundation; either version 2 of the License, or
11  *  (at your option) any later version.
12  *
13  *  This program is distributed in the hope that it will be useful,
14  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
15  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
16  *  GNU General Public License for more details.
17  *
18  *  You should have received a copy of the GNU General Public License
19  *  along with this program; if not, write to the Free Software
20  *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
21  */
22
23 #include <vector>
24 #include <stdexcept>
25
26 #include "mul.h"
27 #include "add.h"
28 #include "power.h"
29 #include "matrix.h"
30 #include "archive.h"
31 #include "utils.h"
32
33 namespace GiNaC {
34
35 GINAC_IMPLEMENT_REGISTERED_CLASS(mul, expairseq)
36
37 //////////
38 // default ctor, dtor, copy ctor, assignment operator and helpers
39 //////////
40
41 mul::mul()
42 {
43         tinfo_key = TINFO_mul;
44 }
45
46 DEFAULT_COPY(mul)
47 DEFAULT_DESTROY(mul)
48
49 //////////
50 // other ctors
51 //////////
52
53 // public
54
55 mul::mul(const ex & lh, const ex & rh)
56 {
57         tinfo_key = TINFO_mul;
58         overall_coeff = _ex1;
59         construct_from_2_ex(lh,rh);
60         GINAC_ASSERT(is_canonical());
61 }
62
63 mul::mul(const exvector & v)
64 {
65         tinfo_key = TINFO_mul;
66         overall_coeff = _ex1;
67         construct_from_exvector(v);
68         GINAC_ASSERT(is_canonical());
69 }
70
71 mul::mul(const epvector & v)
72 {
73         tinfo_key = TINFO_mul;
74         overall_coeff = _ex1;
75         construct_from_epvector(v);
76         GINAC_ASSERT(is_canonical());
77 }
78
79 mul::mul(const epvector & v, const ex & oc)
80 {
81         tinfo_key = TINFO_mul;
82         overall_coeff = oc;
83         construct_from_epvector(v);
84         GINAC_ASSERT(is_canonical());
85 }
86
87 mul::mul(epvector * vp, const ex & oc)
88 {
89         tinfo_key = TINFO_mul;
90         GINAC_ASSERT(vp!=0);
91         overall_coeff = oc;
92         construct_from_epvector(*vp);
93         delete vp;
94         GINAC_ASSERT(is_canonical());
95 }
96
97 mul::mul(const ex & lh, const ex & mh, const ex & rh)
98 {
99         tinfo_key = TINFO_mul;
100         exvector factors;
101         factors.reserve(3);
102         factors.push_back(lh);
103         factors.push_back(mh);
104         factors.push_back(rh);
105         overall_coeff = _ex1;
106         construct_from_exvector(factors);
107         GINAC_ASSERT(is_canonical());
108 }
109
110 //////////
111 // archiving
112 //////////
113
114 DEFAULT_ARCHIVING(mul)
115
116 //////////
117 // functions overriding virtual functions from base classes
118 //////////
119
120 // public
121
122 void mul::print(const print_context & c, unsigned level) const
123 {
124         if (is_a<print_tree>(c)) {
125
126                 inherited::print(c, level);
127
128         } else if (is_a<print_csrc>(c)) {
129
130                 if (precedence() <= level)
131                         c.s << "(";
132
133                 if (!overall_coeff.is_equal(_ex1)) {
134                         overall_coeff.print(c, precedence());
135                         c.s << "*";
136                 }
137
138                 // Print arguments, separated by "*" or "/"
139                 epvector::const_iterator it = seq.begin(), itend = seq.end();
140                 while (it != itend) {
141
142                         // If the first argument is a negative integer power, it gets printed as "1.0/<expr>"
143                         if (it == seq.begin() && ex_to<numeric>(it->coeff).is_integer() && it->coeff.compare(_num0) < 0) {
144                                 if (is_a<print_csrc_cl_N>(c))
145                                         c.s << "recip(";
146                                 else
147                                         c.s << "1.0/";
148                         }
149
150                         // If the exponent is 1 or -1, it is left out
151                         if (it->coeff.compare(_ex1) == 0 || it->coeff.compare(_num_1) == 0)
152                                 it->rest.print(c, precedence());
153                         else {
154                                 // Outer parens around ex needed for broken gcc-2.95 parser:
155                                 (ex(power(it->rest, abs(ex_to<numeric>(it->coeff))))).print(c, level);
156                         }
157
158                         // Separator is "/" for negative integer powers, "*" otherwise
159                         ++it;
160                         if (it != itend) {
161                                 if (ex_to<numeric>(it->coeff).is_integer() && it->coeff.compare(_num0) < 0)
162                                         c.s << "/";
163                                 else
164                                         c.s << "*";
165                         }
166                 }
167
168                 if (precedence() <= level)
169                         c.s << ")";
170
171         } else {
172
173                 if (precedence() <= level) {
174                         if (is_a<print_latex>(c))
175                                 c.s << "{(";
176                         else
177                                 c.s << "(";
178                 }
179
180                 bool first = true;
181
182                 // First print the overall numeric coefficient
183                 numeric coeff = ex_to<numeric>(overall_coeff);
184                 if (coeff.csgn() == -1)
185                         c.s << '-';
186                 if (!coeff.is_equal(_num1) &&
187                         !coeff.is_equal(_num_1)) {
188                         if (coeff.is_rational()) {
189                                 if (coeff.is_negative())
190                                         (-coeff).print(c);
191                                 else
192                                         coeff.print(c);
193                         } else {
194                                 if (coeff.csgn() == -1)
195                                         (-coeff).print(c, precedence());
196                                 else
197                                         coeff.print(c, precedence());
198                         }
199                         if (is_a<print_latex>(c))
200                                 c.s << ' ';
201                         else
202                                 c.s << '*';
203                 }
204
205                 // Then proceed with the remaining factors
206                 epvector::const_iterator it = seq.begin(), itend = seq.end();
207                 while (it != itend) {
208                         if (!first) {
209                                 if (is_a<print_latex>(c))
210                                         c.s << ' ';
211                                 else
212                                         c.s << '*';
213                         } else {
214                                 first = false;
215                         }
216                         recombine_pair_to_ex(*it).print(c, precedence());
217                         ++it;
218                 }
219
220                 if (precedence() <= level) {
221                         if (is_a<print_latex>(c))
222                                 c.s << ")}";
223                         else
224                                 c.s << ")";
225                 }
226         }
227 }
228
229 bool mul::info(unsigned inf) const
230 {
231         switch (inf) {
232                 case info_flags::polynomial:
233                 case info_flags::integer_polynomial:
234                 case info_flags::cinteger_polynomial:
235                 case info_flags::rational_polynomial:
236                 case info_flags::crational_polynomial:
237                 case info_flags::rational_function: {
238                         epvector::const_iterator i = seq.begin(), end = seq.end();
239                         while (i != end) {
240                                 if (!(recombine_pair_to_ex(*i).info(inf)))
241                                         return false;
242                                 ++i;
243                         }
244                         return overall_coeff.info(inf);
245                 }
246                 case info_flags::algebraic: {
247                         epvector::const_iterator i = seq.begin(), end = seq.end();
248                         while (i != end) {
249                                 if ((recombine_pair_to_ex(*i).info(inf)))
250                                         return true;
251                                 ++i;
252                         }
253                         return false;
254                 }
255         }
256         return inherited::info(inf);
257 }
258
259 int mul::degree(const ex & s) const
260 {
261         // Sum up degrees of factors
262         int deg_sum = 0;
263         epvector::const_iterator i = seq.begin(), end = seq.end();
264         while (i != end) {
265                 if (ex_to<numeric>(i->coeff).is_integer())
266                         deg_sum += i->rest.degree(s) * ex_to<numeric>(i->coeff).to_int();
267                 ++i;
268         }
269         return deg_sum;
270 }
271
272 int mul::ldegree(const ex & s) const
273 {
274         // Sum up degrees of factors
275         int deg_sum = 0;
276         epvector::const_iterator i = seq.begin(), end = seq.end();
277         while (i != end) {
278                 if (ex_to<numeric>(i->coeff).is_integer())
279                         deg_sum += i->rest.ldegree(s) * ex_to<numeric>(i->coeff).to_int();
280                 ++i;
281         }
282         return deg_sum;
283 }
284
285 ex mul::coeff(const ex & s, int n) const
286 {
287         exvector coeffseq;
288         coeffseq.reserve(seq.size()+1);
289         
290         if (n==0) {
291                 // product of individual coeffs
292                 // if a non-zero power of s is found, the resulting product will be 0
293                 epvector::const_iterator i = seq.begin(), end = seq.end();
294                 while (i != end) {
295                         coeffseq.push_back(recombine_pair_to_ex(*i).coeff(s,n));
296                         ++i;
297                 }
298                 coeffseq.push_back(overall_coeff);
299                 return (new mul(coeffseq))->setflag(status_flags::dynallocated);
300         }
301         
302         epvector::const_iterator i = seq.begin(), end = seq.end();
303         bool coeff_found = false;
304         while (i != end) {
305                 ex t = recombine_pair_to_ex(*i);
306                 ex c = t.coeff(s, n);
307                 if (!c.is_zero()) {
308                         coeffseq.push_back(c);
309                         coeff_found = 1;
310                 } else {
311                         coeffseq.push_back(t);
312                 }
313                 ++i;
314         }
315         if (coeff_found) {
316                 coeffseq.push_back(overall_coeff);
317                 return (new mul(coeffseq))->setflag(status_flags::dynallocated);
318         }
319         
320         return _ex0;
321 }
322
323 /** Perform automatic term rewriting rules in this class.  In the following
324  *  x, x1, x2,... stand for a symbolic variables of type ex and c, c1, c2...
325  *  stand for such expressions that contain a plain number.
326  *  - *(...,x;0) -> 0
327  *  - *(+(x1,x2,...);c) -> *(+(*(x1,c),*(x2,c),...))
328  *  - *(x;1) -> x
329  *  - *(;c) -> c
330  *
331  *  @param level cut-off in recursive evaluation */
332 ex mul::eval(int level) const
333 {
334         epvector *evaled_seqp = evalchildren(level);
335         if (evaled_seqp) {
336                 // do more evaluation later
337                 return (new mul(evaled_seqp,overall_coeff))->
338                            setflag(status_flags::dynallocated);
339         }
340         
341 #ifdef DO_GINAC_ASSERT
342         epvector::const_iterator i = seq.begin(), end = seq.end();
343         while (i != end) {
344                 GINAC_ASSERT((!is_exactly_a<mul>(i->rest)) ||
345                              (!(ex_to<numeric>(i->coeff).is_integer())));
346                 GINAC_ASSERT(!(i->is_canonical_numeric()));
347                 if (is_ex_exactly_of_type(recombine_pair_to_ex(*i), numeric))
348                     print(print_tree(std::cerr));
349                 GINAC_ASSERT(!is_exactly_a<numeric>(recombine_pair_to_ex(*i)));
350                 /* for paranoia */
351                 expair p = split_ex_to_pair(recombine_pair_to_ex(*i));
352                 GINAC_ASSERT(p.rest.is_equal(i->rest));
353                 GINAC_ASSERT(p.coeff.is_equal(i->coeff));
354                 /* end paranoia */
355                 ++i;
356         }
357 #endif // def DO_GINAC_ASSERT
358         
359         if (flags & status_flags::evaluated) {
360                 GINAC_ASSERT(seq.size()>0);
361                 GINAC_ASSERT(seq.size()>1 || !overall_coeff.is_equal(_ex1));
362                 return *this;
363         }
364         
365         int seq_size = seq.size();
366         if (overall_coeff.is_zero()) {
367                 // *(...,x;0) -> 0
368                 return _ex0;
369         } else if (seq_size==0) {
370                 // *(;c) -> c
371                 return overall_coeff;
372         } else if (seq_size==1 && overall_coeff.is_equal(_ex1)) {
373                 // *(x;1) -> x
374                 return recombine_pair_to_ex(*(seq.begin()));
375         } else if ((seq_size==1) &&
376                    is_ex_exactly_of_type((*seq.begin()).rest,add) &&
377                    ex_to<numeric>((*seq.begin()).coeff).is_equal(_num1)) {
378                 // *(+(x,y,...);c) -> +(*(x,c),*(y,c),...) (c numeric(), no powers of +())
379                 const add & addref = ex_to<add>((*seq.begin()).rest);
380                 epvector *distrseq = new epvector();
381                 distrseq->reserve(addref.seq.size());
382                 epvector::const_iterator i = addref.seq.begin(), end = addref.seq.end();
383                 while (i != end) {
384                         distrseq->push_back(addref.combine_pair_with_coeff_to_pair(*i, overall_coeff));
385                         ++i;
386                 }
387                 return (new add(distrseq,
388                                 ex_to<numeric>(addref.overall_coeff).
389                                 mul_dyn(ex_to<numeric>(overall_coeff))))
390                       ->setflag(status_flags::dynallocated | status_flags::evaluated);
391         }
392         return this->hold();
393 }
394
395 ex mul::evalf(int level) const
396 {
397         if (level==1)
398                 return mul(seq,overall_coeff);
399         
400         if (level==-max_recursion_level)
401                 throw(std::runtime_error("max recursion level reached"));
402         
403         epvector *s = new epvector();
404         s->reserve(seq.size());
405
406         --level;
407         epvector::const_iterator i = seq.begin(), end = seq.end();
408         while (i != end) {
409                 s->push_back(combine_ex_with_coeff_to_pair(i->rest.evalf(level),
410                                                            i->coeff));
411                 ++i;
412         }
413         return mul(s, overall_coeff.evalf(level));
414 }
415
416 ex mul::evalm(void) const
417 {
418         // numeric*matrix
419         if (seq.size() == 1 && seq[0].coeff.is_equal(_ex1)
420          && is_ex_of_type(seq[0].rest, matrix))
421                 return ex_to<matrix>(seq[0].rest).mul(ex_to<numeric>(overall_coeff));
422
423         // Evaluate children first, look whether there are any matrices at all
424         // (there can be either no matrices or one matrix; if there were more
425         // than one matrix, it would be a non-commutative product)
426         epvector *s = new epvector;
427         s->reserve(seq.size());
428
429         bool have_matrix = false;
430         epvector::iterator the_matrix;
431
432         epvector::const_iterator i = seq.begin(), end = seq.end();
433         while (i != end) {
434                 const ex &m = recombine_pair_to_ex(*i).evalm();
435                 s->push_back(split_ex_to_pair(m));
436                 if (is_ex_of_type(m, matrix)) {
437                         have_matrix = true;
438                         the_matrix = s->end() - 1;
439                 }
440                 ++i;
441         }
442
443         if (have_matrix) {
444
445                 // The product contained a matrix. We will multiply all other factors
446                 // into that matrix.
447                 matrix m = ex_to<matrix>(the_matrix->rest);
448                 s->erase(the_matrix);
449                 ex scalar = (new mul(s, overall_coeff))->setflag(status_flags::dynallocated);
450                 return m.mul_scalar(scalar);
451
452         } else
453                 return (new mul(s, overall_coeff))->setflag(status_flags::dynallocated);
454 }
455
456 ex mul::simplify_ncmul(const exvector & v) const
457 {
458         if (seq.empty())
459                 return inherited::simplify_ncmul(v);
460
461         // Find first noncommutative element and call its simplify_ncmul()
462         epvector::const_iterator i = seq.begin(), end = seq.end();
463         while (i != end) {
464                 if (i->rest.return_type() == return_types::noncommutative)
465                         return i->rest.simplify_ncmul(v);
466                 ++i;
467         }
468         return inherited::simplify_ncmul(v);
469 }
470
471 // protected
472
473 /** Implementation of ex::diff() for a product.  It applies the product rule.
474  *  @see ex::diff */
475 ex mul::derivative(const symbol & s) const
476 {
477         unsigned num = seq.size();
478         exvector addseq;
479         addseq.reserve(num);
480         
481         // D(a*b*c) = D(a)*b*c + a*D(b)*c + a*b*D(c)
482         epvector mulseq = seq;
483         epvector::const_iterator i = seq.begin(), end = seq.end();
484         epvector::iterator i2 = mulseq.begin();
485         while (i != end) {
486                 expair ep = split_ex_to_pair(power(i->rest, i->coeff - _ex1) *
487                                              i->rest.diff(s));
488                 ep.swap(*i2);
489                 addseq.push_back((new mul(mulseq, overall_coeff * i->coeff))->setflag(status_flags::dynallocated));
490                 ep.swap(*i2);
491                 ++i; ++i2;
492         }
493         return (new add(addseq))->setflag(status_flags::dynallocated);
494 }
495
496 int mul::compare_same_type(const basic & other) const
497 {
498         return inherited::compare_same_type(other);
499 }
500
501 bool mul::is_equal_same_type(const basic & other) const
502 {
503         return inherited::is_equal_same_type(other);
504 }
505
506 unsigned mul::return_type(void) const
507 {
508         if (seq.empty()) {
509                 // mul without factors: should not happen, but commutes
510                 return return_types::commutative;
511         }
512         
513         bool all_commutative = true;
514         epvector::const_iterator noncommutative_element; // point to first found nc element
515         
516         epvector::const_iterator i = seq.begin(), end = seq.end();
517         while (i != end) {
518                 unsigned rt = i->rest.return_type();
519                 if (rt == return_types::noncommutative_composite)
520                         return rt; // one ncc -> mul also ncc
521                 if ((rt == return_types::noncommutative) && (all_commutative)) {
522                         // first nc element found, remember position
523                         noncommutative_element = i;
524                         all_commutative = false;
525                 }
526                 if ((rt == return_types::noncommutative) && (!all_commutative)) {
527                         // another nc element found, compare type_infos
528                         if (noncommutative_element->rest.return_type_tinfo() != i->rest.return_type_tinfo()) {
529                                 // diffent types -> mul is ncc
530                                 return return_types::noncommutative_composite;
531                         }
532                 }
533                 ++i;
534         }
535         // all factors checked
536         return all_commutative ? return_types::commutative : return_types::noncommutative;
537 }
538    
539 unsigned mul::return_type_tinfo(void) const
540 {
541         if (seq.empty())
542                 return tinfo_key;  // mul without factors: should not happen
543         
544         // return type_info of first noncommutative element
545         epvector::const_iterator i = seq.begin(), end = seq.end();
546         while (i != end) {
547                 if (i->rest.return_type() == return_types::noncommutative)
548                         return i->rest.return_type_tinfo();
549                 ++i;
550         }
551         // no noncommutative element found, should not happen
552         return tinfo_key;
553 }
554
555 ex mul::thisexpairseq(const epvector & v, const ex & oc) const
556 {
557         return (new mul(v, oc))->setflag(status_flags::dynallocated);
558 }
559
560 ex mul::thisexpairseq(epvector * vp, const ex & oc) const
561 {
562         return (new mul(vp, oc))->setflag(status_flags::dynallocated);
563 }
564
565 expair mul::split_ex_to_pair(const ex & e) const
566 {
567         if (is_ex_exactly_of_type(e,power)) {
568                 const power & powerref = ex_to<power>(e);
569                 if (is_ex_exactly_of_type(powerref.exponent,numeric))
570                         return expair(powerref.basis,powerref.exponent);
571         }
572         return expair(e,_ex1);
573 }
574         
575 expair mul::combine_ex_with_coeff_to_pair(const ex & e,
576                                           const ex & c) const
577 {
578         // to avoid duplication of power simplification rules,
579         // we create a temporary power object
580         // otherwise it would be hard to correctly simplify
581         // expression like (4^(1/3))^(3/2)
582         if (are_ex_trivially_equal(c,_ex1))
583                 return split_ex_to_pair(e);
584         
585         return split_ex_to_pair(power(e,c));
586 }
587         
588 expair mul::combine_pair_with_coeff_to_pair(const expair & p,
589                                             const ex & c) const
590 {
591         // to avoid duplication of power simplification rules,
592         // we create a temporary power object
593         // otherwise it would be hard to correctly simplify
594         // expression like (4^(1/3))^(3/2)
595         if (are_ex_trivially_equal(c,_ex1))
596                 return p;
597         
598         return split_ex_to_pair(power(recombine_pair_to_ex(p),c));
599 }
600         
601 ex mul::recombine_pair_to_ex(const expair & p) const
602 {
603         if (ex_to<numeric>(p.coeff).is_equal(_num1)) 
604                 return p.rest;
605         else
606                 return power(p.rest,p.coeff);
607 }
608
609 bool mul::expair_needs_further_processing(epp it)
610 {
611         if (is_ex_exactly_of_type((*it).rest,mul) &&
612                 ex_to<numeric>((*it).coeff).is_integer()) {
613                 // combined pair is product with integer power -> expand it
614                 *it = split_ex_to_pair(recombine_pair_to_ex(*it));
615                 return true;
616         }
617         if (is_ex_exactly_of_type((*it).rest,numeric)) {
618                 expair ep=split_ex_to_pair(recombine_pair_to_ex(*it));
619                 if (!ep.is_equal(*it)) {
620                         // combined pair is a numeric power which can be simplified
621                         *it = ep;
622                         return true;
623                 }
624                 if (ex_to<numeric>((*it).coeff).is_equal(_num1)) {
625                         // combined pair has coeff 1 and must be moved to the end
626                         return true;
627                 }
628         }
629         return false;
630 }       
631
632 ex mul::default_overall_coeff(void) const
633 {
634         return _ex1;
635 }
636
637 void mul::combine_overall_coeff(const ex & c)
638 {
639         GINAC_ASSERT(is_exactly_a<numeric>(overall_coeff));
640         GINAC_ASSERT(is_exactly_a<numeric>(c));
641         overall_coeff = ex_to<numeric>(overall_coeff).mul_dyn(ex_to<numeric>(c));
642 }
643
644 void mul::combine_overall_coeff(const ex & c1, const ex & c2)
645 {
646         GINAC_ASSERT(is_exactly_a<numeric>(overall_coeff));
647         GINAC_ASSERT(is_exactly_a<numeric>(c1));
648         GINAC_ASSERT(is_exactly_a<numeric>(c2));
649         overall_coeff = ex_to<numeric>(overall_coeff).mul_dyn(ex_to<numeric>(c1).power(ex_to<numeric>(c2)));
650 }
651
652 bool mul::can_make_flat(const expair & p) const
653 {
654         GINAC_ASSERT(is_exactly_a<numeric>(p.coeff));
655         // this assertion will probably fail somewhere
656         // it would require a more careful make_flat, obeying the power laws
657         // probably should return true only if p.coeff is integer
658         return ex_to<numeric>(p.coeff).is_equal(_num1);
659 }
660
661 ex mul::expand(unsigned options) const
662 {
663         // First, expand the children
664         epvector * expanded_seqp = expandchildren(options);
665         const epvector & expanded_seq = (expanded_seqp == NULL) ? seq : *expanded_seqp;
666
667         // Now, look for all the factors that are sums and multiply each one out
668         // with the next one that is found while collecting the factors which are
669         // not sums
670         int number_of_adds = 0;
671         ex last_expanded = _ex1;
672         epvector non_adds;
673         non_adds.reserve(expanded_seq.size());
674         epvector::const_iterator cit = expanded_seq.begin(), last = expanded_seq.end();
675         while (cit != last) {
676                 if (is_ex_exactly_of_type(cit->rest, add) &&
677                         (cit->coeff.is_equal(_ex1))) {
678                         ++number_of_adds;
679                         if (is_ex_exactly_of_type(last_expanded, add)) {
680                                 const add & add1 = ex_to<add>(last_expanded);
681                                 const add & add2 = ex_to<add>(cit->rest);
682                                 int n1 = add1.nops();
683                                 int n2 = add2.nops();
684                                 exvector distrseq;
685                                 distrseq.reserve(n1*n2);
686                                 for (int i1=0; i1<n1; ++i1) {
687                                         for (int i2=0; i2<n2; ++i2) {
688                                                 distrseq.push_back(add1.op(i1) * add2.op(i2));
689                                         }
690                                 }
691                                 last_expanded = (new add(distrseq))->
692                                                  setflag(status_flags::dynallocated | (options == 0 ? status_flags::expanded : 0));
693                         } else {
694                                 non_adds.push_back(split_ex_to_pair(last_expanded));
695                                 last_expanded = cit->rest;
696                         }
697                 } else {
698                         non_adds.push_back(*cit);
699                 }
700                 ++cit;
701         }
702         if (expanded_seqp)
703                 delete expanded_seqp;
704         
705         // Now the only remaining thing to do is to multiply the factors which
706         // were not sums into the "last_expanded" sum
707         if (is_ex_exactly_of_type(last_expanded, add)) {
708                 const add & finaladd = ex_to<add>(last_expanded);
709                 exvector distrseq;
710                 int n = finaladd.nops();
711                 distrseq.reserve(n);
712                 for (int i=0; i<n; ++i) {
713                         epvector factors = non_adds;
714                         factors.push_back(split_ex_to_pair(finaladd.op(i)));
715                         distrseq.push_back((new mul(factors, overall_coeff))->
716                                             setflag(status_flags::dynallocated | (options == 0 ? status_flags::expanded : 0)));
717                 }
718                 return ((new add(distrseq))->
719                         setflag(status_flags::dynallocated | (options == 0 ? status_flags::expanded : 0)));
720         }
721         non_adds.push_back(split_ex_to_pair(last_expanded));
722         return (new mul(non_adds, overall_coeff))->
723                 setflag(status_flags::dynallocated | (options == 0 ? status_flags::expanded : 0));
724 }
725
726   
727 //////////
728 // new virtual functions which can be overridden by derived classes
729 //////////
730
731 // none
732
733 //////////
734 // non-virtual functions in this class
735 //////////
736
737
738 /** Member-wise expand the expairs representing this sequence.  This must be
739  *  overridden from expairseq::expandchildren() and done iteratively in order
740  *  to allow for early cancallations and thus safe memory.
741  *
742  *  @see mul::expand()
743  *  @return pointer to epvector containing expanded representation or zero
744  *  pointer, if sequence is unchanged. */
745 epvector * mul::expandchildren(unsigned options) const
746 {
747         epvector::const_iterator last = seq.end();
748         epvector::const_iterator cit = seq.begin();
749         while (cit!=last) {
750                 const ex & factor = recombine_pair_to_ex(*cit);
751                 const ex & expanded_factor = factor.expand(options);
752                 if (!are_ex_trivially_equal(factor,expanded_factor)) {
753                         
754                         // something changed, copy seq, eval and return it
755                         epvector *s = new epvector;
756                         s->reserve(seq.size());
757                         
758                         // copy parts of seq which are known not to have changed
759                         epvector::const_iterator cit2 = seq.begin();
760                         while (cit2!=cit) {
761                                 s->push_back(*cit2);
762                                 ++cit2;
763                         }
764                         // copy first changed element
765                         s->push_back(split_ex_to_pair(expanded_factor));
766                         ++cit2;
767                         // copy rest
768                         while (cit2!=last) {
769                                 s->push_back(split_ex_to_pair(recombine_pair_to_ex(*cit2).expand(options)));
770                                 ++cit2;
771                         }
772                         return s;
773                 }
774                 ++cit;
775         }
776         
777         return 0; // nothing has changed
778 }
779
780 } // namespace GiNaC