- Inserted some more std:: to make it compile under GCC2.96.
[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                 deg_sum+=(*cit).rest.degree(s) * ex_to_numeric((*cit).coeff).to_int();
314         }
315         return deg_sum;
316 }
317
318 int mul::ldegree(const symbol & s) const
319 {
320         int deg_sum = 0;
321         for (epvector::const_iterator cit=seq.begin(); cit!=seq.end(); ++cit) {
322                 deg_sum+=(*cit).rest.ldegree(s) * ex_to_numeric((*cit).coeff).to_int();
323         }
324         return deg_sum;
325 }
326
327 ex mul::coeff(const symbol & s, int n) const
328 {
329         exvector coeffseq;
330         coeffseq.reserve(seq.size()+1);
331         
332         if (n==0) {
333                 // product of individual coeffs
334                 // if a non-zero power of s is found, the resulting product will be 0
335                 epvector::const_iterator it=seq.begin();
336                 while (it!=seq.end()) {
337                         coeffseq.push_back(recombine_pair_to_ex(*it).coeff(s,n));
338                         ++it;
339                 }
340                 coeffseq.push_back(overall_coeff);
341                 return (new mul(coeffseq))->setflag(status_flags::dynallocated);
342         }
343                  
344         epvector::const_iterator it=seq.begin();
345         bool coeff_found=0;
346         while (it!=seq.end()) {
347                 ex t=recombine_pair_to_ex(*it);
348                 ex c=t.coeff(s,n);
349                 if (!c.is_zero()) {
350                         coeffseq.push_back(c);
351                         coeff_found=1;
352                 } else {
353                         coeffseq.push_back(t);
354                 }
355                 ++it;
356         }
357         if (coeff_found) {
358                 coeffseq.push_back(overall_coeff);
359                 return (new mul(coeffseq))->setflag(status_flags::dynallocated);
360         }
361         
362         return _ex0();
363 }
364
365 ex mul::eval(int level) const
366 {
367         // simplifications  *(...,x;0) -> 0
368         //                  *(+(x,y,...);c) -> *(+(*(x,c),*(y,c),...)) (c numeric())
369         //                  *(x;1) -> x
370         //                  *(;c) -> c
371
372         debugmsg("mul eval",LOGLEVEL_MEMBER_FUNCTION);
373
374         epvector * evaled_seqp=evalchildren(level);
375         if (evaled_seqp!=0) {
376                 // do more evaluation later
377                 return (new mul(evaled_seqp,overall_coeff))->
378                                    setflag(status_flags::dynallocated);
379         }
380
381 #ifdef DO_GINAC_ASSERT
382         for (epvector::const_iterator cit=seq.begin(); cit!=seq.end(); ++cit) {
383                 GINAC_ASSERT((!is_ex_exactly_of_type((*cit).rest,mul))||
384                            (!(ex_to_numeric((*cit).coeff).is_integer())));
385                 GINAC_ASSERT(!((*cit).is_numeric_with_coeff_1()));
386                 if (is_ex_exactly_of_type(recombine_pair_to_ex(*cit),numeric)) {
387                         printtree(cerr,0);
388                 }
389                 GINAC_ASSERT(!is_ex_exactly_of_type(recombine_pair_to_ex(*cit),numeric));
390                 /* for paranoia */
391                 expair p=split_ex_to_pair(recombine_pair_to_ex(*cit));
392                 GINAC_ASSERT(p.rest.is_equal((*cit).rest));
393                 GINAC_ASSERT(p.coeff.is_equal((*cit).coeff));
394                 /* end paranoia */
395         }
396 #endif // def DO_GINAC_ASSERT
397
398         if (flags & status_flags::evaluated) {
399                 GINAC_ASSERT(seq.size()>0);
400                 GINAC_ASSERT((seq.size()>1)||!overall_coeff.is_equal(_ex1()));
401                 return *this;
402         }
403
404         int seq_size=seq.size();
405         if (overall_coeff.is_equal(_ex0())) {
406                 // *(...,x;0) -> 0
407                 return _ex0();
408         } else if (seq_size==0) {
409                 // *(;c) -> c
410                 return overall_coeff;
411         } else if ((seq_size==1)&&overall_coeff.is_equal(_ex1())) {
412                 // *(x;1) -> x
413                 return recombine_pair_to_ex(*(seq.begin()));
414         } else if ((seq_size==1) &&
415                    is_ex_exactly_of_type((*seq.begin()).rest,add) &&
416                    ex_to_numeric((*seq.begin()).coeff).is_equal(_num1())) {
417                 // *(+(x,y,...);c) -> +(*(x,c),*(y,c),...) (c numeric(), no powers of +())
418                 const add & addref=ex_to_add((*seq.begin()).rest);
419                 epvector distrseq;
420                 distrseq.reserve(addref.seq.size());
421                 for (epvector::const_iterator cit=addref.seq.begin(); cit!=addref.seq.end(); ++cit) {
422                         distrseq.push_back(addref.combine_pair_with_coeff_to_pair(*cit, overall_coeff));
423                 }
424                 return (new add(distrseq,
425                                 ex_to_numeric(addref.overall_coeff).
426                                 mul_dyn(ex_to_numeric(overall_coeff))))
427                       ->setflag(status_flags::dynallocated | status_flags::evaluated);
428         }
429         return this->hold();
430 }
431
432 ex mul::evalf(int level) const
433 {
434         if (level==1)
435                 return mul(seq,overall_coeff);
436         
437         if (level==-max_recursion_level)
438                 throw(std::runtime_error("max recursion level reached"));
439         
440         epvector s;
441         s.reserve(seq.size());
442         
443         --level;
444         for (epvector::const_iterator it=seq.begin(); it!=seq.end(); ++it) {
445                 s.push_back(combine_ex_with_coeff_to_pair((*it).rest.evalf(level),
446                                                           (*it).coeff));
447         }
448         return mul(s,overall_coeff.evalf(level));
449 }
450
451 exvector mul::get_indices(void) const
452 {
453         // return union of indices of factors
454         exvector iv;
455         for (epvector::const_iterator cit=seq.begin(); cit!=seq.end(); ++cit) {
456                 exvector subiv=(*cit).rest.get_indices();
457                 iv.reserve(iv.size()+subiv.size());
458                 for (exvector::const_iterator cit2=subiv.begin(); cit2!=subiv.end(); ++cit2) {
459                         iv.push_back(*cit2);
460                 }
461         }
462         return iv;
463 }
464
465 ex mul::simplify_ncmul(const exvector & v) const
466 {
467         throw(std::logic_error("mul::simplify_ncmul() should never have been called!"));
468 }
469
470 // protected
471
472 /** Implementation of ex::diff() for a product.  It applies the product rule.
473  *  @see ex::diff */
474 ex mul::derivative(const symbol & s) const
475 {
476         exvector addseq;
477         addseq.reserve(seq.size());
478         
479         // D(a*b*c) = D(a)*b*c + a*D(b)*c + a*b*D(c)
480         for (unsigned i=0; i!=seq.size(); ++i) {
481                 epvector mulseq = seq;
482                 mulseq[i] = split_ex_to_pair(power(seq[i].rest,seq[i].coeff - _ex1()) *
483                                              seq[i].rest.diff(s));
484                 addseq.push_back((new mul(mulseq,overall_coeff*seq[i].coeff))->setflag(status_flags::dynallocated));
485         }
486         return (new add(addseq))->setflag(status_flags::dynallocated);
487 }
488
489 int mul::compare_same_type(const basic & other) const
490 {
491         return inherited::compare_same_type(other);
492 }
493
494 bool mul::is_equal_same_type(const basic & other) const
495 {
496         return inherited::is_equal_same_type(other);
497 }
498
499 unsigned mul::return_type(void) const
500 {
501         if (seq.size()==0) {
502                 // mul without factors: should not happen, but commutes
503                 return return_types::commutative;
504         }
505
506         bool all_commutative = 1;
507         unsigned rt;
508         epvector::const_iterator cit_noncommutative_element; // point to first found nc element
509
510         for (epvector::const_iterator cit=seq.begin(); cit!=seq.end(); ++cit) {
511                 rt=(*cit).rest.return_type();
512                 if (rt==return_types::noncommutative_composite) return rt; // one ncc -> mul also ncc
513                 if ((rt==return_types::noncommutative)&&(all_commutative)) {
514                         // first nc element found, remember position
515                         cit_noncommutative_element = cit;
516                         all_commutative = 0;
517                 }
518                 if ((rt==return_types::noncommutative)&&(!all_commutative)) {
519                         // another nc element found, compare type_infos
520                         if ((*cit_noncommutative_element).rest.return_type_tinfo()!=(*cit).rest.return_type_tinfo()) {
521                                 // diffent types -> mul is ncc
522                                 return return_types::noncommutative_composite;
523                         }
524                 }
525         }
526         // all factors checked
527         return all_commutative ? return_types::commutative : return_types::noncommutative;
528 }
529    
530 unsigned mul::return_type_tinfo(void) const
531 {
532         if (seq.size()==0) {
533                 // mul without factors: should not happen
534                 return tinfo_key;
535         }
536         // return type_info of first noncommutative element
537         for (epvector::const_iterator cit=seq.begin(); cit!=seq.end(); ++cit) {
538                 if ((*cit).rest.return_type()==return_types::noncommutative) {
539                         return (*cit).rest.return_type_tinfo();
540                 }
541         }
542         // no noncommutative element found, should not happen
543         return tinfo_key;
544 }
545
546 ex mul::thisexpairseq(const epvector & v, const ex & oc) const
547 {
548         return (new mul(v,oc))->setflag(status_flags::dynallocated);
549 }
550
551 ex mul::thisexpairseq(epvector * vp, const ex & oc) const
552 {
553         return (new mul(vp,oc))->setflag(status_flags::dynallocated);
554 }
555
556 expair mul::split_ex_to_pair(const ex & e) const
557 {
558         if (is_ex_exactly_of_type(e,power)) {
559                 const power & powerref=ex_to_power(e);
560                 if (is_ex_exactly_of_type(powerref.exponent,numeric)) {
561                         return expair(powerref.basis,powerref.exponent);
562                 }
563         }
564         return expair(e,_ex1());
565 }
566         
567 expair mul::combine_ex_with_coeff_to_pair(const ex & e,
568                                           const ex & c) const
569 {
570         // to avoid duplication of power simplification rules,
571         // we create a temporary power object
572         // otherwise it would be hard to correctly simplify
573         // expression like (4^(1/3))^(3/2)
574         if (are_ex_trivially_equal(c,_ex1()))
575                 return split_ex_to_pair(e);
576         
577         return split_ex_to_pair(power(e,c));
578 }
579         
580 expair mul::combine_pair_with_coeff_to_pair(const expair & p,
581                                             const ex & c) const
582 {
583         // to avoid duplication of power simplification rules,
584         // we create a temporary power object
585         // otherwise it would be hard to correctly simplify
586         // expression like (4^(1/3))^(3/2)
587         if (are_ex_trivially_equal(c,_ex1()))
588                 return p;
589         
590         return split_ex_to_pair(power(recombine_pair_to_ex(p),c));
591 }
592         
593 ex mul::recombine_pair_to_ex(const expair & p) const
594 {
595         if (ex_to_numeric(p.coeff).is_equal(_num1())) 
596                 return p.rest;
597         else
598                 return power(p.rest,p.coeff);
599 }
600
601 bool mul::expair_needs_further_processing(epp it)
602 {
603         if (is_ex_exactly_of_type((*it).rest,mul) &&
604                 ex_to_numeric((*it).coeff).is_integer()) {
605                 // combined pair is product with integer power -> expand it
606                 *it=split_ex_to_pair(recombine_pair_to_ex(*it));
607                 return true;
608         }
609         if (is_ex_exactly_of_type((*it).rest,numeric)) {
610                 expair ep=split_ex_to_pair(recombine_pair_to_ex(*it));
611                 if (!ep.is_equal(*it)) {
612                         // combined pair is a numeric power which can be simplified
613                         *it=ep;
614                         return true;
615                 }
616                 if (ex_to_numeric((*it).coeff).is_equal(_num1())) {
617                         // combined pair has coeff 1 and must be moved to the end
618                         return true;
619                 }
620         }
621         return false;
622 }       
623
624 ex mul::default_overall_coeff(void) const
625 {
626         return _ex1();
627 }
628
629 void mul::combine_overall_coeff(const ex & c)
630 {
631         GINAC_ASSERT(is_ex_exactly_of_type(overall_coeff,numeric));
632         GINAC_ASSERT(is_ex_exactly_of_type(c,numeric));
633         overall_coeff = ex_to_numeric(overall_coeff).mul_dyn(ex_to_numeric(c));
634 }
635
636 void mul::combine_overall_coeff(const ex & c1, const ex & c2)
637 {
638         GINAC_ASSERT(is_ex_exactly_of_type(overall_coeff,numeric));
639         GINAC_ASSERT(is_ex_exactly_of_type(c1,numeric));
640         GINAC_ASSERT(is_ex_exactly_of_type(c2,numeric));
641         overall_coeff = ex_to_numeric(overall_coeff).mul_dyn(ex_to_numeric(c1).power(ex_to_numeric(c2)));
642 }
643
644 bool mul::can_make_flat(const expair & p) const
645 {
646         GINAC_ASSERT(is_ex_exactly_of_type(p.coeff,numeric));
647         // this assertion will probably fail somewhere
648         // it would require a more careful make_flat, obeying the power laws
649         // probably should return true only if p.coeff is integer
650         return ex_to_numeric(p.coeff).is_equal(_num1());
651 }
652
653 ex mul::expand(unsigned options) const
654 {
655         if (flags & status_flags::expanded)
656                 return *this;
657         
658         exvector sub_expanded_seq;
659         intvector positions_of_adds;
660         intvector number_of_add_operands;
661         
662         epvector * expanded_seqp = expandchildren(options);
663         
664         const epvector & expanded_seq = expanded_seqp==0 ? seq : *expanded_seqp;
665         
666         positions_of_adds.resize(expanded_seq.size());
667         number_of_add_operands.resize(expanded_seq.size());
668         
669         int number_of_adds = 0;
670         int number_of_expanded_terms = 1;
671         
672         unsigned current_position = 0;
673         epvector::const_iterator last = expanded_seq.end();
674         for (epvector::const_iterator cit = expanded_seq.begin(); cit!=last; ++cit) {
675                 if (is_ex_exactly_of_type((*cit).rest,add) &&
676                         ((*cit).coeff.is_equal(_ex1()))) {
677                         positions_of_adds[number_of_adds] = current_position;
678                         const add & expanded_addref = ex_to_add((*cit).rest);
679                         unsigned addref_nops = expanded_addref.nops();
680                         number_of_add_operands[number_of_adds] = addref_nops;
681                         number_of_expanded_terms *= addref_nops;
682                         ++number_of_adds;
683                 }
684                 ++current_position;
685         }
686         
687         if (number_of_adds==0) {
688                 if (expanded_seqp==0)
689                         return this->setflag(status_flags::expanded);
690                 else
691                         return ((new mul(expanded_seqp,overall_coeff))->
692                                setflag(status_flags::dynallocated | status_flags::expanded));
693         }
694         
695         exvector distrseq;
696         distrseq.reserve(number_of_expanded_terms);
697         
698         intvector k;
699         k.resize(number_of_adds, 0);
700         
701         for (;;) {
702                 epvector term;
703                 term = expanded_seq;
704                 for (int l=0; l<number_of_adds; ++l) {
705                         const add & addref = ex_to_add(expanded_seq[positions_of_adds[l]].rest);
706                         GINAC_ASSERT(term[positions_of_adds[l]].coeff.compare(_ex1())==0);
707                         term[positions_of_adds[l]]=split_ex_to_pair(addref.op(k[l]));
708                 }
709                 distrseq.push_back((new mul(term,overall_coeff))->
710                                     setflag(status_flags::dynallocated | status_flags::expanded));
711                 
712                 // increment k[]
713                 int l = number_of_adds-1;
714                 while ((l>=0) && ((++k[l])>=number_of_add_operands[l])) {
715                         k[l] = 0;    
716                         --l;
717                 }
718                 if (l < 0) break;
719         }
720         
721         if (expanded_seqp!=0)
722                 delete expanded_seqp;
723         
724         return (new add(distrseq))->setflag(status_flags::dynallocated |
725                                                                                 status_flags::expanded);
726 }
727
728 //////////
729 // new virtual functions which can be overridden by derived classes
730 //////////
731
732 // none
733
734 //////////
735 // non-virtual functions in this class
736 //////////
737
738 epvector * mul::expandchildren(unsigned options) const
739 {
740         epvector::const_iterator last = seq.end();
741         epvector::const_iterator cit = seq.begin();
742         while (cit!=last) {
743                 const ex & factor = recombine_pair_to_ex(*cit);
744                 const ex & expanded_factor = factor.expand(options);
745                 if (!are_ex_trivially_equal(factor,expanded_factor)) {
746                         
747                         // something changed, copy seq, eval and return it
748                         epvector *s=new epvector;
749                         s->reserve(seq.size());
750                         
751                         // copy parts of seq which are known not to have changed
752                         epvector::const_iterator cit2 = seq.begin();
753                         while (cit2!=cit) {
754                                 s->push_back(*cit2);
755                                 ++cit2;
756                         }
757                         // copy first changed element
758                         s->push_back(split_ex_to_pair(expanded_factor));
759                         ++cit2;
760                         // copy rest
761                         while (cit2!=last) {
762                                 s->push_back(split_ex_to_pair(recombine_pair_to_ex(*cit2).expand(options)));
763                                 ++cit2;
764                         }
765                         return s;
766                 }
767                 ++cit;
768         }
769         
770         return 0; // nothing has changed
771 }
772    
773 //////////
774 // static member variables
775 //////////
776
777 // protected
778
779 unsigned mul::precedence = 50;
780
781
782 //////////
783 // global constants
784 //////////
785
786 const mul some_mul;
787 const std::type_info & typeid_mul = typeid(some_mul);
788
789 #ifndef NO_NAMESPACE_GINAC
790 } // namespace GiNaC
791 #endif // ndef NO_NAMESPACE_GINAC