- non-integer powers are treated as 0 by (l)degree() and coeff()
[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-2000 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 "archive.h"
30 #include "debugmsg.h"
31 #include "utils.h"
32
33 #ifndef NO_NAMESPACE_GINAC
34 namespace GiNaC {
35 #endif // ndef NO_NAMESPACE_GINAC
36
37 GINAC_IMPLEMENT_REGISTERED_CLASS(mul, expairseq)
38
39 //////////
40 // default constructor, destructor, copy constructor assignment operator and helpers
41 //////////
42
43 // public
44
45 mul::mul()
46 {
47         debugmsg("mul default constructor",LOGLEVEL_CONSTRUCT);
48         tinfo_key = TINFO_mul;
49 }
50
51 mul::~mul()
52 {
53         debugmsg("mul destructor",LOGLEVEL_DESTRUCT);
54         destroy(false);
55 }
56
57 mul::mul(const mul & other)
58 {
59         debugmsg("mul copy constructor",LOGLEVEL_CONSTRUCT);
60         copy(other);
61 }
62
63 const mul & mul::operator=(const mul & other)
64 {
65         debugmsg("mul operator=",LOGLEVEL_ASSIGNMENT);
66         if (this != &other) {
67                 destroy(true);
68                 copy(other);
69         }
70         return *this;
71 }
72
73 // protected
74
75 void mul::copy(const mul & other)
76 {
77         inherited::copy(other);
78 }
79
80 void mul::destroy(bool call_parent)
81 {
82         if (call_parent) inherited::destroy(call_parent);
83 }
84
85 //////////
86 // other constructors
87 //////////
88
89 // public
90
91 mul::mul(const ex & lh, const ex & rh)
92 {
93         debugmsg("mul constructor from ex,ex",LOGLEVEL_CONSTRUCT);
94         tinfo_key = TINFO_mul;
95         overall_coeff = _ex1();
96         construct_from_2_ex(lh,rh);
97         GINAC_ASSERT(is_canonical());
98 }
99
100 mul::mul(const exvector & v)
101 {
102         debugmsg("mul constructor from exvector",LOGLEVEL_CONSTRUCT);
103         tinfo_key = TINFO_mul;
104         overall_coeff = _ex1();
105         construct_from_exvector(v);
106         GINAC_ASSERT(is_canonical());
107 }
108
109 mul::mul(const epvector & v)
110 {
111         debugmsg("mul constructor from epvector",LOGLEVEL_CONSTRUCT);
112         tinfo_key = TINFO_mul;
113         overall_coeff = _ex1();
114         construct_from_epvector(v);
115         GINAC_ASSERT(is_canonical());
116 }
117
118 mul::mul(const epvector & v, const ex & oc)
119 {
120         debugmsg("mul constructor from epvector,ex",LOGLEVEL_CONSTRUCT);
121         tinfo_key = TINFO_mul;
122         overall_coeff = oc;
123         construct_from_epvector(v);
124         GINAC_ASSERT(is_canonical());
125 }
126
127 mul::mul(epvector * vp, const ex & oc)
128 {
129         debugmsg("mul constructor from epvector *,ex",LOGLEVEL_CONSTRUCT);
130         tinfo_key = TINFO_mul;
131         GINAC_ASSERT(vp!=0);
132         overall_coeff = oc;
133         construct_from_epvector(*vp);
134         delete vp;
135         GINAC_ASSERT(is_canonical());
136 }
137
138 mul::mul(const ex & lh, const ex & mh, const ex & rh)
139 {
140         debugmsg("mul constructor from ex,ex,ex",LOGLEVEL_CONSTRUCT);
141         tinfo_key = TINFO_mul;
142         exvector factors;
143         factors.reserve(3);
144         factors.push_back(lh);
145         factors.push_back(mh);
146         factors.push_back(rh);
147         overall_coeff = _ex1();
148         construct_from_exvector(factors);
149         GINAC_ASSERT(is_canonical());
150 }
151
152 //////////
153 // archiving
154 //////////
155
156 /** Construct object from archive_node. */
157 mul::mul(const archive_node &n, const lst &sym_lst) : inherited(n, sym_lst)
158 {
159         debugmsg("mul constructor from archive_node", LOGLEVEL_CONSTRUCT);
160 }
161
162 /** Unarchive the object. */
163 ex mul::unarchive(const archive_node &n, const lst &sym_lst)
164 {
165         return (new mul(n, sym_lst))->setflag(status_flags::dynallocated);
166 }
167
168 /** Archive the object. */
169 void mul::archive(archive_node &n) const
170 {
171         inherited::archive(n);
172 }
173
174 //////////
175 // functions overriding virtual functions from bases classes
176 //////////
177
178 // public
179
180 basic * mul::duplicate() const
181 {
182         debugmsg("mul duplicate",LOGLEVEL_ASSIGNMENT);
183         return new mul(*this);
184 }
185
186 void mul::print(std::ostream & os, unsigned upper_precedence) const
187 {
188         debugmsg("mul print",LOGLEVEL_PRINT);
189         if (precedence<=upper_precedence) os << "(";
190         bool first=true;
191         // first print the overall numeric coefficient:
192         numeric coeff = ex_to_numeric(overall_coeff);
193         if (coeff.csgn()==-1) os << '-';
194         if (!coeff.is_equal(_num1()) &&
195                 !coeff.is_equal(_num_1())) {
196                 if (coeff.is_rational()) {
197                         if (coeff.is_negative())
198                                 os << -coeff;
199                         else
200                                 os << coeff;
201                 } else {
202                         if (coeff.csgn()==-1)
203                                 (-coeff).print(os, precedence);
204                         else
205                                 coeff.print(os, precedence);
206                 }
207                 os << '*';
208         }
209         // then proceed with the remaining factors:
210         for (epvector::const_iterator cit=seq.begin(); cit!=seq.end(); ++cit) {
211                 if (!first) {
212                         os << '*';
213                 } else {
214                         first=false;
215                 }
216                 recombine_pair_to_ex(*cit).print(os,precedence);
217         }
218         if (precedence<=upper_precedence) os << ")";
219 }
220
221 void mul::printraw(std::ostream & os) const
222 {
223         debugmsg("mul printraw",LOGLEVEL_PRINT);
224
225         os << "*(";
226         for (epvector::const_iterator it=seq.begin(); it!=seq.end(); ++it) {
227                 os << "(";
228                 (*it).rest.bp->printraw(os);
229                 os << ",";
230                 (*it).coeff.bp->printraw(os);
231                 os << "),";
232         }
233         os << ",hash=" << hashvalue << ",flags=" << flags;
234         os << ")";
235 }
236
237 void mul::printcsrc(std::ostream & os, unsigned type, unsigned upper_precedence) const
238 {
239         debugmsg("mul print csrc", LOGLEVEL_PRINT);
240         if (precedence <= upper_precedence)
241                 os << "(";
242
243         if (!overall_coeff.is_equal(_ex1())) {
244                 overall_coeff.bp->printcsrc(os,type,precedence);
245                 os << "*";
246         }
247         
248         // Print arguments, separated by "*" or "/"
249         epvector::const_iterator it = seq.begin();
250         epvector::const_iterator itend = seq.end();
251         while (it != itend) {
252
253                 // If the first argument is a negative integer power, it gets printed as "1.0/<expr>"
254                 if (it == seq.begin() && ex_to_numeric(it->coeff).is_integer() && it->coeff.compare(_num0()) < 0) {
255                         if (type == csrc_types::ctype_cl_N)
256                                 os << "recip(";
257                         else
258                                 os << "1.0/";
259                 }
260
261                 // If the exponent is 1 or -1, it is left out
262                 if (it->coeff.compare(_ex1()) == 0 || it->coeff.compare(_num_1()) == 0)
263                         it->rest.bp->printcsrc(os, type, precedence);
264                 else
265                         // outer parens around ex needed for broken gcc-2.95 parser:
266                         (ex(power(it->rest, abs(ex_to_numeric(it->coeff))))).bp->printcsrc(os, type, upper_precedence);
267
268                 // Separator is "/" for negative integer powers, "*" otherwise
269                 ++it;
270                 if (it != itend) {
271                         if (ex_to_numeric(it->coeff).is_integer() && it->coeff.compare(_num0()) < 0)
272                                 os << "/";
273                         else
274                                 os << "*";
275                 }
276         }
277         if (precedence <= upper_precedence)
278                 os << ")";
279 }
280
281 bool mul::info(unsigned inf) const
282 {
283         switch (inf) {
284                 case info_flags::polynomial:
285                 case info_flags::integer_polynomial:
286                 case info_flags::cinteger_polynomial:
287                 case info_flags::rational_polynomial:
288                 case info_flags::crational_polynomial:
289                 case info_flags::rational_function: {
290                         for (epvector::const_iterator i=seq.begin(); i!=seq.end(); ++i) {
291                                 if (!(recombine_pair_to_ex(*i).info(inf)))
292                                         return false;
293                         }
294                         return overall_coeff.info(inf);
295                 }
296                 case info_flags::algebraic: {
297                         for (epvector::const_iterator i=seq.begin(); i!=seq.end(); ++i) {
298                                 if ((recombine_pair_to_ex(*i).info(inf)))
299                                         return true;
300                         }
301                         return false;
302                 }
303         }
304         return inherited::info(inf);
305 }
306
307 typedef std::vector<int> intvector;
308
309 int mul::degree(const symbol & s) const
310 {
311         int deg_sum = 0;
312         for (epvector::const_iterator cit=seq.begin(); cit!=seq.end(); ++cit) {
313                 if (ex_to_numeric(cit->coeff).is_integer())
314                         deg_sum+=cit->rest.degree(s) * ex_to_numeric(cit->coeff).to_int();
315         }
316         return deg_sum;
317 }
318
319 int mul::ldegree(const symbol & s) const
320 {
321         int deg_sum = 0;
322         for (epvector::const_iterator cit=seq.begin(); cit!=seq.end(); ++cit) {
323                 if (ex_to_numeric(cit->coeff).is_integer())
324                         deg_sum+=cit->rest.ldegree(s) * ex_to_numeric(cit->coeff).to_int();
325         }
326         return deg_sum;
327 }
328
329 ex mul::coeff(const symbol & s, int n) const
330 {
331         exvector coeffseq;
332         coeffseq.reserve(seq.size()+1);
333         
334         if (n==0) {
335                 // product of individual coeffs
336                 // if a non-zero power of s is found, the resulting product will be 0
337                 epvector::const_iterator it=seq.begin();
338                 while (it!=seq.end()) {
339                         coeffseq.push_back(recombine_pair_to_ex(*it).coeff(s,n));
340                         ++it;
341                 }
342                 coeffseq.push_back(overall_coeff);
343                 return (new mul(coeffseq))->setflag(status_flags::dynallocated);
344         }
345                  
346         epvector::const_iterator it=seq.begin();
347         bool coeff_found=0;
348         while (it!=seq.end()) {
349                 ex t=recombine_pair_to_ex(*it);
350                 ex c=t.coeff(s,n);
351                 if (!c.is_zero()) {
352                         coeffseq.push_back(c);
353                         coeff_found=1;
354                 } else {
355                         coeffseq.push_back(t);
356                 }
357                 ++it;
358         }
359         if (coeff_found) {
360                 coeffseq.push_back(overall_coeff);
361                 return (new mul(coeffseq))->setflag(status_flags::dynallocated);
362         }
363         
364         return _ex0();
365 }
366
367 ex mul::eval(int level) const
368 {
369         // simplifications  *(...,x;0) -> 0
370         //                  *(+(x,y,...);c) -> *(+(*(x,c),*(y,c),...)) (c numeric())
371         //                  *(x;1) -> x
372         //                  *(;c) -> c
373
374         debugmsg("mul eval",LOGLEVEL_MEMBER_FUNCTION);
375
376         epvector * evaled_seqp=evalchildren(level);
377         if (evaled_seqp!=0) {
378                 // do more evaluation later
379                 return (new mul(evaled_seqp,overall_coeff))->
380                                    setflag(status_flags::dynallocated);
381         }
382
383 #ifdef DO_GINAC_ASSERT
384         for (epvector::const_iterator cit=seq.begin(); cit!=seq.end(); ++cit) {
385                 GINAC_ASSERT((!is_ex_exactly_of_type((*cit).rest,mul))||
386                            (!(ex_to_numeric((*cit).coeff).is_integer())));
387                 GINAC_ASSERT(!((*cit).is_numeric_with_coeff_1()));
388                 if (is_ex_exactly_of_type(recombine_pair_to_ex(*cit),numeric)) {
389                         printtree(cerr,0);
390                 }
391                 GINAC_ASSERT(!is_ex_exactly_of_type(recombine_pair_to_ex(*cit),numeric));
392                 /* for paranoia */
393                 expair p=split_ex_to_pair(recombine_pair_to_ex(*cit));
394                 GINAC_ASSERT(p.rest.is_equal((*cit).rest));
395                 GINAC_ASSERT(p.coeff.is_equal((*cit).coeff));
396                 /* end paranoia */
397         }
398 #endif // def DO_GINAC_ASSERT
399
400         if (flags & status_flags::evaluated) {
401                 GINAC_ASSERT(seq.size()>0);
402                 GINAC_ASSERT((seq.size()>1)||!overall_coeff.is_equal(_ex1()));
403                 return *this;
404         }
405
406         int seq_size=seq.size();
407         if (overall_coeff.is_equal(_ex0())) {
408                 // *(...,x;0) -> 0
409                 return _ex0();
410         } else if (seq_size==0) {
411                 // *(;c) -> c
412                 return overall_coeff;
413         } else if ((seq_size==1)&&overall_coeff.is_equal(_ex1())) {
414                 // *(x;1) -> x
415                 return recombine_pair_to_ex(*(seq.begin()));
416         } else if ((seq_size==1) &&
417                    is_ex_exactly_of_type((*seq.begin()).rest,add) &&
418                    ex_to_numeric((*seq.begin()).coeff).is_equal(_num1())) {
419                 // *(+(x,y,...);c) -> +(*(x,c),*(y,c),...) (c numeric(), no powers of +())
420                 const add & addref=ex_to_add((*seq.begin()).rest);
421                 epvector distrseq;
422                 distrseq.reserve(addref.seq.size());
423                 for (epvector::const_iterator cit=addref.seq.begin(); cit!=addref.seq.end(); ++cit) {
424                         distrseq.push_back(addref.combine_pair_with_coeff_to_pair(*cit, overall_coeff));
425                 }
426                 return (new add(distrseq,
427                                 ex_to_numeric(addref.overall_coeff).
428                                 mul_dyn(ex_to_numeric(overall_coeff))))
429                       ->setflag(status_flags::dynallocated | status_flags::evaluated);
430         }
431         return this->hold();
432 }
433
434 ex mul::evalf(int level) const
435 {
436         if (level==1)
437                 return mul(seq,overall_coeff);
438         
439         if (level==-max_recursion_level)
440                 throw(std::runtime_error("max recursion level reached"));
441         
442         epvector s;
443         s.reserve(seq.size());
444         
445         --level;
446         for (epvector::const_iterator it=seq.begin(); it!=seq.end(); ++it) {
447                 s.push_back(combine_ex_with_coeff_to_pair((*it).rest.evalf(level),
448                                                           (*it).coeff));
449         }
450         return mul(s,overall_coeff.evalf(level));
451 }
452
453 exvector mul::get_indices(void) const
454 {
455         // return union of indices of factors
456         exvector iv;
457         for (epvector::const_iterator cit=seq.begin(); cit!=seq.end(); ++cit) {
458                 exvector subiv=(*cit).rest.get_indices();
459                 iv.reserve(iv.size()+subiv.size());
460                 for (exvector::const_iterator cit2=subiv.begin(); cit2!=subiv.end(); ++cit2) {
461                         iv.push_back(*cit2);
462                 }
463         }
464         return iv;
465 }
466
467 ex mul::simplify_ncmul(const exvector & v) const
468 {
469         throw(std::logic_error("mul::simplify_ncmul() should never have been called!"));
470 }
471
472 // protected
473
474 /** Implementation of ex::diff() for a product.  It applies the product rule.
475  *  @see ex::diff */
476 ex mul::derivative(const symbol & s) const
477 {
478         exvector addseq;
479         addseq.reserve(seq.size());
480         
481         // D(a*b*c) = D(a)*b*c + a*D(b)*c + a*b*D(c)
482         for (unsigned i=0; i!=seq.size(); ++i) {
483                 epvector mulseq = seq;
484                 mulseq[i] = split_ex_to_pair(power(seq[i].rest,seq[i].coeff - _ex1()) *
485                                              seq[i].rest.diff(s));
486                 addseq.push_back((new mul(mulseq,overall_coeff*seq[i].coeff))->setflag(status_flags::dynallocated));
487         }
488         return (new add(addseq))->setflag(status_flags::dynallocated);
489 }
490
491 int mul::compare_same_type(const basic & other) const
492 {
493         return inherited::compare_same_type(other);
494 }
495
496 bool mul::is_equal_same_type(const basic & other) const
497 {
498         return inherited::is_equal_same_type(other);
499 }
500
501 unsigned mul::return_type(void) const
502 {
503         if (seq.size()==0) {
504                 // mul without factors: should not happen, but commutes
505                 return return_types::commutative;
506         }
507
508         bool all_commutative = 1;
509         unsigned rt;
510         epvector::const_iterator cit_noncommutative_element; // point to first found nc element
511
512         for (epvector::const_iterator cit=seq.begin(); cit!=seq.end(); ++cit) {
513                 rt=(*cit).rest.return_type();
514                 if (rt==return_types::noncommutative_composite) return rt; // one ncc -> mul also ncc
515                 if ((rt==return_types::noncommutative)&&(all_commutative)) {
516                         // first nc element found, remember position
517                         cit_noncommutative_element = cit;
518                         all_commutative = 0;
519                 }
520                 if ((rt==return_types::noncommutative)&&(!all_commutative)) {
521                         // another nc element found, compare type_infos
522                         if ((*cit_noncommutative_element).rest.return_type_tinfo()!=(*cit).rest.return_type_tinfo()) {
523                                 // diffent types -> mul is ncc
524                                 return return_types::noncommutative_composite;
525                         }
526                 }
527         }
528         // all factors checked
529         return all_commutative ? return_types::commutative : return_types::noncommutative;
530 }
531    
532 unsigned mul::return_type_tinfo(void) const
533 {
534         if (seq.size()==0) {
535                 // mul without factors: should not happen
536                 return tinfo_key;
537         }
538         // return type_info of first noncommutative element
539         for (epvector::const_iterator cit=seq.begin(); cit!=seq.end(); ++cit) {
540                 if ((*cit).rest.return_type()==return_types::noncommutative) {
541                         return (*cit).rest.return_type_tinfo();
542                 }
543         }
544         // no noncommutative element found, should not happen
545         return tinfo_key;
546 }
547
548 ex mul::thisexpairseq(const epvector & v, const ex & oc) const
549 {
550         return (new mul(v,oc))->setflag(status_flags::dynallocated);
551 }
552
553 ex mul::thisexpairseq(epvector * vp, const ex & oc) const
554 {
555         return (new mul(vp,oc))->setflag(status_flags::dynallocated);
556 }
557
558 expair mul::split_ex_to_pair(const ex & e) const
559 {
560         if (is_ex_exactly_of_type(e,power)) {
561                 const power & powerref=ex_to_power(e);
562                 if (is_ex_exactly_of_type(powerref.exponent,numeric)) {
563                         return expair(powerref.basis,powerref.exponent);
564                 }
565         }
566         return expair(e,_ex1());
567 }
568         
569 expair mul::combine_ex_with_coeff_to_pair(const ex & e,
570                                           const ex & c) const
571 {
572         // to avoid duplication of power simplification rules,
573         // we create a temporary power object
574         // otherwise it would be hard to correctly simplify
575         // expression like (4^(1/3))^(3/2)
576         if (are_ex_trivially_equal(c,_ex1()))
577                 return split_ex_to_pair(e);
578         
579         return split_ex_to_pair(power(e,c));
580 }
581         
582 expair mul::combine_pair_with_coeff_to_pair(const expair & p,
583                                             const ex & c) const
584 {
585         // to avoid duplication of power simplification rules,
586         // we create a temporary power object
587         // otherwise it would be hard to correctly simplify
588         // expression like (4^(1/3))^(3/2)
589         if (are_ex_trivially_equal(c,_ex1()))
590                 return p;
591         
592         return split_ex_to_pair(power(recombine_pair_to_ex(p),c));
593 }
594         
595 ex mul::recombine_pair_to_ex(const expair & p) const
596 {
597         if (ex_to_numeric(p.coeff).is_equal(_num1())) 
598                 return p.rest;
599         else
600                 return power(p.rest,p.coeff);
601 }
602
603 bool mul::expair_needs_further_processing(epp it)
604 {
605         if (is_ex_exactly_of_type((*it).rest,mul) &&
606                 ex_to_numeric((*it).coeff).is_integer()) {
607                 // combined pair is product with integer power -> expand it
608                 *it=split_ex_to_pair(recombine_pair_to_ex(*it));
609                 return true;
610         }
611         if (is_ex_exactly_of_type((*it).rest,numeric)) {
612                 expair ep=split_ex_to_pair(recombine_pair_to_ex(*it));
613                 if (!ep.is_equal(*it)) {
614                         // combined pair is a numeric power which can be simplified
615                         *it=ep;
616                         return true;
617                 }
618                 if (ex_to_numeric((*it).coeff).is_equal(_num1())) {
619                         // combined pair has coeff 1 and must be moved to the end
620                         return true;
621                 }
622         }
623         return false;
624 }       
625
626 ex mul::default_overall_coeff(void) const
627 {
628         return _ex1();
629 }
630
631 void mul::combine_overall_coeff(const ex & c)
632 {
633         GINAC_ASSERT(is_ex_exactly_of_type(overall_coeff,numeric));
634         GINAC_ASSERT(is_ex_exactly_of_type(c,numeric));
635         overall_coeff = ex_to_numeric(overall_coeff).mul_dyn(ex_to_numeric(c));
636 }
637
638 void mul::combine_overall_coeff(const ex & c1, const ex & c2)
639 {
640         GINAC_ASSERT(is_ex_exactly_of_type(overall_coeff,numeric));
641         GINAC_ASSERT(is_ex_exactly_of_type(c1,numeric));
642         GINAC_ASSERT(is_ex_exactly_of_type(c2,numeric));
643         overall_coeff = ex_to_numeric(overall_coeff).mul_dyn(ex_to_numeric(c1).power(ex_to_numeric(c2)));
644 }
645
646 bool mul::can_make_flat(const expair & p) const
647 {
648         GINAC_ASSERT(is_ex_exactly_of_type(p.coeff,numeric));
649         // this assertion will probably fail somewhere
650         // it would require a more careful make_flat, obeying the power laws
651         // probably should return true only if p.coeff is integer
652         return ex_to_numeric(p.coeff).is_equal(_num1());
653 }
654
655 ex mul::expand(unsigned options) const
656 {
657         if (flags & status_flags::expanded)
658                 return *this;
659         
660         exvector sub_expanded_seq;
661         intvector positions_of_adds;
662         intvector number_of_add_operands;
663         
664         epvector * expanded_seqp = expandchildren(options);
665         
666         const epvector & expanded_seq = expanded_seqp==0 ? seq : *expanded_seqp;
667         
668         positions_of_adds.resize(expanded_seq.size());
669         number_of_add_operands.resize(expanded_seq.size());
670         
671         int number_of_adds = 0;
672         int number_of_expanded_terms = 1;
673         
674         unsigned current_position = 0;
675         epvector::const_iterator last = expanded_seq.end();
676         for (epvector::const_iterator cit = expanded_seq.begin(); cit!=last; ++cit) {
677                 if (is_ex_exactly_of_type((*cit).rest,add) &&
678                         ((*cit).coeff.is_equal(_ex1()))) {
679                         positions_of_adds[number_of_adds] = current_position;
680                         const add & expanded_addref = ex_to_add((*cit).rest);
681                         unsigned addref_nops = expanded_addref.nops();
682                         number_of_add_operands[number_of_adds] = addref_nops;
683                         number_of_expanded_terms *= addref_nops;
684                         ++number_of_adds;
685                 }
686                 ++current_position;
687         }
688         
689         if (number_of_adds==0) {
690                 if (expanded_seqp==0)
691                         return this->setflag(status_flags::expanded);
692                 else
693                         return ((new mul(expanded_seqp,overall_coeff))->
694                                setflag(status_flags::dynallocated | status_flags::expanded));
695         }
696         
697         exvector distrseq;
698         distrseq.reserve(number_of_expanded_terms);
699         
700         intvector k;
701         k.resize(number_of_adds, 0);
702         
703         for (;;) {
704                 epvector term;
705                 term = expanded_seq;
706                 for (int l=0; l<number_of_adds; ++l) {
707                         const add & addref = ex_to_add(expanded_seq[positions_of_adds[l]].rest);
708                         GINAC_ASSERT(term[positions_of_adds[l]].coeff.compare(_ex1())==0);
709                         term[positions_of_adds[l]]=split_ex_to_pair(addref.op(k[l]));
710                 }
711                 distrseq.push_back((new mul(term,overall_coeff))->
712                                     setflag(status_flags::dynallocated | status_flags::expanded));
713                 
714                 // increment k[]
715                 int l = number_of_adds-1;
716                 while ((l>=0) && ((++k[l])>=number_of_add_operands[l])) {
717                         k[l] = 0;    
718                         --l;
719                 }
720                 if (l < 0) break;
721         }
722         
723         if (expanded_seqp!=0)
724                 delete expanded_seqp;
725         
726         return (new add(distrseq))->setflag(status_flags::dynallocated |
727                                                                                 status_flags::expanded);
728 }
729
730 //////////
731 // new virtual functions which can be overridden by derived classes
732 //////////
733
734 // none
735
736 //////////
737 // non-virtual functions in this class
738 //////////
739
740 epvector * mul::expandchildren(unsigned options) const
741 {
742         epvector::const_iterator last = seq.end();
743         epvector::const_iterator cit = seq.begin();
744         while (cit!=last) {
745                 const ex & factor = recombine_pair_to_ex(*cit);
746                 const ex & expanded_factor = factor.expand(options);
747                 if (!are_ex_trivially_equal(factor,expanded_factor)) {
748                         
749                         // something changed, copy seq, eval and return it
750                         epvector *s=new epvector;
751                         s->reserve(seq.size());
752                         
753                         // copy parts of seq which are known not to have changed
754                         epvector::const_iterator cit2 = seq.begin();
755                         while (cit2!=cit) {
756                                 s->push_back(*cit2);
757                                 ++cit2;
758                         }
759                         // copy first changed element
760                         s->push_back(split_ex_to_pair(expanded_factor));
761                         ++cit2;
762                         // copy rest
763                         while (cit2!=last) {
764                                 s->push_back(split_ex_to_pair(recombine_pair_to_ex(*cit2).expand(options)));
765                                 ++cit2;
766                         }
767                         return s;
768                 }
769                 ++cit;
770         }
771         
772         return 0; // nothing has changed
773 }
774    
775 //////////
776 // static member variables
777 //////////
778
779 // protected
780
781 unsigned mul::precedence = 50;
782
783
784 //////////
785 // global constants
786 //////////
787
788 const mul some_mul;
789 const std::type_info & typeid_mul = typeid(some_mul);
790
791 #ifndef NO_NAMESPACE_GINAC
792 } // namespace GiNaC
793 #endif // ndef NO_NAMESPACE_GINAC