5054fb06af5871d5e04bd5ab6a61198dda2817b1
[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 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 "debugmsg.h"
30
31 namespace GiNaC {
32
33 //////////
34 // default constructor, destructor, copy constructor assignment operator and helpers
35 //////////
36
37 // public
38
39 mul::mul()
40 {
41     debugmsg("mul default constructor",LOGLEVEL_CONSTRUCT);
42     tinfo_key = TINFO_mul;
43 }
44
45 mul::~mul()
46 {
47     debugmsg("mul destructor",LOGLEVEL_DESTRUCT);
48     destroy(0);
49 }
50
51 mul::mul(mul const & other)
52 {
53     debugmsg("mul copy constructor",LOGLEVEL_CONSTRUCT);
54     copy(other);
55 }
56
57 mul const & mul::operator=(mul const & other)
58 {
59     debugmsg("mul operator=",LOGLEVEL_ASSIGNMENT);
60     if (this != &other) {
61         destroy(1);
62         copy(other);
63     }
64     return *this;
65 }
66
67 // protected
68
69 void mul::copy(mul const & other)
70 {
71     expairseq::copy(other);
72 }
73
74 void mul::destroy(bool call_parent)
75 {
76     if (call_parent) expairseq::destroy(call_parent);
77 }
78
79 //////////
80 // other constructors
81 //////////
82
83 // public
84
85 mul::mul(ex const & lh, ex const & rh)
86 {
87     debugmsg("mul constructor from ex,ex",LOGLEVEL_CONSTRUCT);
88     tinfo_key = TINFO_mul;
89     overall_coeff=exONE();
90     construct_from_2_ex(lh,rh);
91     ASSERT(is_canonical());
92 }
93
94 mul::mul(exvector const & v)
95 {
96     debugmsg("mul constructor from exvector",LOGLEVEL_CONSTRUCT);
97     tinfo_key = TINFO_mul;
98     overall_coeff=exONE();
99     construct_from_exvector(v);
100     ASSERT(is_canonical());
101 }
102
103 /*
104 mul::mul(epvector const & v, bool do_not_canonicalize)
105 {
106     debugmsg("mul constructor from epvector,bool",LOGLEVEL_CONSTRUCT);
107     tinfo_key = TINFO_mul;
108     if (do_not_canonicalize) {
109         seq=v;
110 #ifdef EXPAIRSEQ_USE_HASHTAB
111         combine_same_terms(); // to build hashtab
112 #endif // def EXPAIRSEQ_USE_HASHTAB
113     } else {
114         construct_from_epvector(v);
115     }
116     ASSERT(is_canonical());
117 }
118 */
119
120 mul::mul(epvector const & v)
121 {
122     debugmsg("mul constructor from epvector",LOGLEVEL_CONSTRUCT);
123     tinfo_key = TINFO_mul;
124     overall_coeff=exONE();
125     construct_from_epvector(v);
126     ASSERT(is_canonical());
127 }
128
129 mul::mul(epvector const & v, ex const & oc)
130 {
131     debugmsg("mul constructor from epvector,ex",LOGLEVEL_CONSTRUCT);
132     tinfo_key = TINFO_mul;
133     overall_coeff=oc;
134     construct_from_epvector(v);
135     ASSERT(is_canonical());
136 }
137
138 mul::mul(epvector * vp, ex const & oc)
139 {
140     debugmsg("mul constructor from epvector *,ex",LOGLEVEL_CONSTRUCT);
141     tinfo_key = TINFO_mul;
142     ASSERT(vp!=0);
143     overall_coeff=oc;
144     construct_from_epvector(*vp);
145     delete vp;
146     ASSERT(is_canonical());
147 }
148
149 mul::mul(ex const & lh, ex const & mh, ex const & rh)
150 {
151     debugmsg("mul constructor from ex,ex,ex",LOGLEVEL_CONSTRUCT);
152     tinfo_key = TINFO_mul;
153     exvector factors;
154     factors.reserve(3);
155     factors.push_back(lh);
156     factors.push_back(mh);
157     factors.push_back(rh);
158     overall_coeff=exONE();
159     construct_from_exvector(factors);
160     ASSERT(is_canonical());
161 }
162
163 //////////
164 // functions overriding virtual functions from bases classes
165 //////////
166
167 // public
168
169 basic * mul::duplicate() const
170 {
171     debugmsg("mul duplicate",LOGLEVEL_ASSIGNMENT);
172     return new mul(*this);
173 }
174
175 bool mul::info(unsigned inf) const
176 {
177     // TODO: optimize
178     if (inf==info_flags::polynomial || inf==info_flags::integer_polynomial || inf==info_flags::rational_polynomial || inf==info_flags::rational_function) {
179         for (epvector::const_iterator it=seq.begin(); it!=seq.end(); ++it) {
180             if (!(recombine_pair_to_ex(*it).info(inf)))
181                 return false;
182         }
183         return true;
184     } else {
185         return expairseq::info(inf);
186     }
187 }
188
189 typedef vector<int> intvector;
190
191 int mul::degree(symbol const & s) const
192 {
193     int deg_sum=0;
194     for (epvector::const_iterator cit=seq.begin(); cit!=seq.end(); ++cit) {
195         deg_sum+=(*cit).rest.degree(s) * ex_to_numeric((*cit).coeff).to_int();
196     }
197     return deg_sum;
198 }
199
200 int mul::ldegree(symbol const & s) const
201 {
202     int deg_sum=0;
203     for (epvector::const_iterator cit=seq.begin(); cit!=seq.end(); ++cit) {
204         deg_sum+=(*cit).rest.ldegree(s) * ex_to_numeric((*cit).coeff).to_int();
205     }
206     return deg_sum;
207 }
208
209 ex mul::coeff(symbol const & s, int const n) const
210 {
211     exvector coeffseq;
212     coeffseq.reserve(seq.size()+1);
213     
214     if (n==0) {
215         // product of individual coeffs
216         // if a non-zero power of s is found, the resulting product will be 0
217         epvector::const_iterator it=seq.begin();
218         while (it!=seq.end()) {
219             coeffseq.push_back(recombine_pair_to_ex(*it).coeff(s,n));
220             ++it;
221         }
222         coeffseq.push_back(overall_coeff);
223         return (new mul(coeffseq))->setflag(status_flags::dynallocated);
224     }
225          
226     epvector::const_iterator it=seq.begin();
227     bool coeff_found=0;
228     while (it!=seq.end()) {
229         ex t=recombine_pair_to_ex(*it);
230         ex c=t.coeff(s,n);
231         if (!c.is_zero()) {
232             coeffseq.push_back(c);
233             coeff_found=1;
234         } else {
235             coeffseq.push_back(t);
236         }
237         ++it;
238     }
239     if (coeff_found) {
240         coeffseq.push_back(overall_coeff);
241         return (new mul(coeffseq))->setflag(status_flags::dynallocated);
242     }
243     
244     return exZERO();
245 }
246
247 /*
248 ex mul::eval(int level) const
249 {
250     // simplifications: *(...,x,(c1,1),(c2,1)) -> *(...,x,(c1*c2,1)) (c1, c2 numeric(), move pairs to end first)
251     //                  *(...,x,1) -> *(...,x)
252     //                  *(...,x,0) -> 0
253     //                  *(+(x,y,...),(c,1)) -> *(+(*(x,c),*(y,c),...)) (c numeric())
254     //                  *(x) -> x
255     //                  *() -> 1
256
257     debugmsg("mul eval",LOGLEVEL_MEMBER_FUNCTION);
258
259     if ((level==1)&&(flags & status_flags::evaluated)) {
260 #ifdef DOASSERT
261         for (epvector::const_iterator cit=seq.begin(); cit!=seq.end(); ++cit) {
262             ASSERT((!is_ex_exactly_of_type((*cit).rest,mul))||
263                    (!(ex_to_numeric((*cit).coeff).is_integer())));
264         }
265
266         // test if all numerics were moved to the end and
267         // all numerics with coeff 1 to the very end
268         if (seq.size()!=0) {
269             epvector::const_iterator cit=seq.end();
270             bool all_coeff_1=true;
271             bool all_numeric=true;
272             do {
273                 cit--;
274                 if (is_ex_exactly_of_type((*cit).rest,numeric)) {
275                     ASSERT(all_numeric);
276                     if ((*cit).coeff.is_equal(exONE())) {
277                         ASSERT(all_coeff_1);
278                     } else {
279                         all_coeff_1=false;
280                     }
281                 } else {
282                     all_numeric=false;
283                 }
284             } while (cit!=seq.begin());
285         }
286 #endif // def DOASSERT    
287         return *this;
288     }
289
290     epvector newseq;
291     epvector::iterator it1,it2;
292     bool seq_copied=false;
293
294     epvector * evaled_seqp=evalchildren(level);
295     if (evaled_seqp!=0) {
296         // do more evaluation later
297         return (new mul(evaled_seqp))->setflag(status_flags::dynallocated);
298     }
299
300     // combine pairs with coeff 1 (all numerics should be at end, assert below)
301     if (seq.size()>1) {
302         // count number of pairs with coeff 1
303         unsigned num_coeff_1=0;
304         bool still_numeric=true;
305         epvector::const_iterator cit=seq.end();
306         unsigned first_pos;
307         unsigned second_pos;
308         do {
309             cit--;
310             if (is_ex_exactly_of_type((*cit).rest,numeric)) {
311                 if ((*cit).coeff.is_equal(exONE())) {
312                     num_coeff_1++;
313                 }
314             } else {
315                 still_numeric=false;
316             }
317         } while ((cit!=seq.begin())&&still_numeric);
318         if (num_coeff_1>1) {
319             newseq=seq;
320             
321     }
322     
323     
324 #ifdef DOASSERT
325     for (epvector::const_iterator cit=seq.begin(); cit!=seq.end(); ++cit) {
326         ASSERT((!is_ex_exactly_of_type((*cit).rest,mul))||
327                (!(ex_to_numeric((*cit).coeff).is_integer())));
328     }
329
330     // test if all numerics were moved to the end and
331     // all numerics with coeff 1 to the very end
332     if (seq.size()!=0) {
333         epvector::const_iterator cit=seq.end();
334         bool all_coeff_1=true;
335         bool all_numeric=true;
336         do {
337             cit--;
338             if (is_ex_exactly_of_type((*cit).rest,numeric)) {
339                 ASSERT(all_numeric);
340                 if ((*cit).coeff.is_equal(exONE())) {
341                     ASSERT(all_coeff_1);
342                 } else {
343                     all_coeff_1=false;
344                 }
345             } else {
346                 all_numeric=false;
347             }
348         } while (cit!=seq.begin());
349     }
350 #endif // def DOASSERT
351     
352     if (flags & status_flags::evaluated) {
353         return *this;
354     }
355     
356     expair const & last_expair=*(seq.end()-1);
357     expair const & next_to_last_expair=*(seq.end()-2);
358     int seq_size = seq.size();
359
360     // *(...,x,(c1,1),(c2,1)) -> *(...,x,(c1*c2,1)) (c1, c2 numeric())
361     if ((!seq_copied) && (seq_size>=2) &&
362         is_ex_exactly_of_type(last_expair.rest,numeric) &&
363         ex_to_numeric(last_expair.coeff).is_equal(numONE()) &&
364         is_ex_exactly_of_type(next_to_last_expair.rest,numeric) &&
365         ex_to_numeric(next_to_last_expair.coeff).is_equal(numONE()) ) {
366         newseq=seq;
367         seq_copied=true;
368         it2=newseq.end()-1;
369         it1=it2-1;
370     }
371     while (seq_copied && (newseq.size()>=2) &&
372            is_ex_exactly_of_type((*it1).rest,numeric) &&
373            ex_to_numeric((*it1).coeff).is_equal(numONE()) &&
374            is_ex_exactly_of_type((*it2).rest,numeric) &&
375            ex_to_numeric((*it2).coeff).is_equal(numONE()) ) {
376         *it1=expair(ex_to_numeric((*it1).rest).mul_dyn(ex_to_numeric((*it2).rest)),exONE());
377         newseq.pop_back();
378         it2=newseq.end()-1;
379         it1=it2-1;
380     }
381
382     // *(...,x,1) -> *(...,x)
383     if ((!seq_copied) && (seq_size>=1) &&
384         (is_ex_exactly_of_type(last_expair.rest,numeric)) &&
385         (ex_to_numeric(last_expair.rest).compare(numONE())==0)) {
386         newseq=seq;
387         seq_copied=true;
388         it2=newseq.end()-1;
389     }
390     if (seq_copied && (newseq.size()>=1) &&
391         (is_ex_exactly_of_type((*it2).rest,numeric)) &&
392         (ex_to_numeric((*it2).rest).compare(numONE())==0)) {
393         newseq.pop_back();
394         it2=newseq.end()-1;
395     }
396
397     // *(...,x,0) -> 0
398     if ((!seq_copied) && (seq_size>=1) &&
399         (is_ex_exactly_of_type(last_expair.rest,numeric)) &&
400         (ex_to_numeric(last_expair.rest).is_zero())) {
401         return exZERO();
402     }
403     if (seq_copied && (newseq.size()>=1) &&
404         (is_ex_exactly_of_type((*it2).rest,numeric)) &&
405         (ex_to_numeric((*it2).rest).is_zero())) {
406         return exZERO();
407     }
408
409     // *(+(x,y,...),c) -> +(*(x,c),*(y,c),...) (c numeric(), no powers of +())
410     if ((!seq_copied) && (seq_size==2) &&
411         is_ex_exactly_of_type(next_to_last_expair.rest,add) &&
412         is_ex_exactly_of_type(last_expair.rest,numeric) &&
413         ex_to_numeric(last_expair.coeff).is_equal(numONE()) &&
414         (ex_to_numeric(next_to_last_expair.coeff).compare(numONE())==0)) {
415         add const & addref=ex_to_add(next_to_last_expair.rest);
416         epvector distrseq;
417         distrseq.reserve(addref.seq.size());
418         for (epvector::const_iterator cit=addref.seq.begin(); cit!=addref.seq.end(); ++cit) {
419             distrseq.push_back(addref.combine_pair_with_coeff_to_pair(*cit,
420                                    last_expair.rest));
421         }
422         // special treatment for the last element if it is numeric (to
423         // avoid terms like (2/3)*(3/2)) is no longer necessary, this
424         // is handled in add::combine_pair_with_coeff_to_pair()
425         return (new add(distrseq,1))->setflag(status_flags::dynallocated  |
426                                               status_flags::evaluated );
427     }
428     if (seq_copied && (newseq.size()==2) &&
429         is_ex_exactly_of_type(newseq[0].rest,add) &&
430         is_ex_exactly_of_type(newseq[1].rest,numeric) &&
431         ex_to_numeric(newseq[1].coeff).is_equal(numONE()) &&
432         (ex_to_numeric(newseq[0].coeff).compare(numONE())==0)) {
433         add const & addref=ex_to_add(newseq[0].rest);
434         epvector distrseq;
435         distrseq.reserve(addref.seq.size());
436         for (epvector::const_iterator cit=addref.seq.begin(); cit!=addref.seq.end(); ++cit) {
437             distrseq.push_back(addref.combine_pair_with_coeff_to_pair(*cit,
438                                    newseq[1].rest));
439         }
440         // special treatment for the last element if it is numeric (to
441         // avoid terms like (2/3)*(3/2)) is no longer necessary, this
442         // is handled in add::combine_pair_with_coeff_to_pair()
443         return (new add(distrseq,1))->setflag(status_flags::dynallocated  |
444                                               status_flags::evaluated );
445     }
446     
447     // *() -> 1
448     if ((!seq_copied) && (seq_size==0)) {
449         return exONE();
450     } else if (seq_copied && (newseq.size()==0)) {
451         return exONE();
452     }
453
454     // *(x) -> x
455     if ((!seq_copied) && (seq_size==1)) {
456         return recombine_pair_to_ex(*(seq.begin()));
457     } else if (seq_copied && (newseq.size()==1)) {
458         return recombine_pair_to_ex(*(newseq.begin()));
459     }
460
461     if (!seq_copied) return this->hold();
462
463     return (new mul(newseq,1))->setflag(status_flags::dynallocated  |
464                                         status_flags::evaluated );
465 }
466 */
467
468 ex mul::eval(int level) const
469 {
470     // simplifications  *(...,x;0) -> 0
471     //                  *(+(x,y,...);c) -> *(+(*(x,c),*(y,c),...)) (c numeric())
472     //                  *(x;1) -> x
473     //                  *(;c) -> c
474
475     debugmsg("mul eval",LOGLEVEL_MEMBER_FUNCTION);
476
477     epvector * evaled_seqp=evalchildren(level);
478     if (evaled_seqp!=0) {
479         // do more evaluation later
480         return (new mul(evaled_seqp,overall_coeff))->
481                    setflag(status_flags::dynallocated);
482     }
483
484 #ifdef DOASSERT
485     for (epvector::const_iterator cit=seq.begin(); cit!=seq.end(); ++cit) {
486         ASSERT((!is_ex_exactly_of_type((*cit).rest,mul))||
487                (!(ex_to_numeric((*cit).coeff).is_integer())));
488         ASSERT(!((*cit).is_numeric_with_coeff_1()));
489         if (is_ex_exactly_of_type(recombine_pair_to_ex(*cit),numeric)) {
490             printtree(cerr,0);
491         }
492         ASSERT(!is_ex_exactly_of_type(recombine_pair_to_ex(*cit),numeric));
493         /* for paranoia */
494         expair p=split_ex_to_pair(recombine_pair_to_ex(*cit));
495         ASSERT(p.rest.is_equal((*cit).rest));
496         ASSERT(p.coeff.is_equal((*cit).coeff));
497         /* end paranoia */
498     }
499 #endif // def DOASSERT
500
501     if (flags & status_flags::evaluated) {
502         ASSERT(seq.size()>0);
503         ASSERT((seq.size()>1)||!overall_coeff.is_equal(exONE()));
504         return *this;
505     }
506
507     int seq_size=seq.size();
508     if (overall_coeff.is_equal(exZERO())) {
509         // *(...,x;0) -> 0
510         return exZERO();
511     } else if (seq_size==0) {
512         // *(;c) -> c
513         return overall_coeff;
514     } else if ((seq_size==1)&&overall_coeff.is_equal(exONE())) {
515         // *(x;1) -> x
516         return recombine_pair_to_ex(*(seq.begin()));
517     } else if ((seq_size==1) &&
518                is_ex_exactly_of_type((*seq.begin()).rest,add) &&
519                ex_to_numeric((*seq.begin()).coeff).is_equal(numONE())) {
520         // *(+(x,y,...);c) -> +(*(x,c),*(y,c),...) (c numeric(), no powers of +())
521         add const & addref=ex_to_add((*seq.begin()).rest);
522         epvector distrseq;
523         distrseq.reserve(addref.seq.size());
524         for (epvector::const_iterator cit=addref.seq.begin(); cit!=addref.seq.end(); ++cit) {
525             distrseq.push_back(addref.combine_pair_with_coeff_to_pair(*cit,
526                                    overall_coeff));
527         }
528         return (new add(distrseq,
529                         ex_to_numeric(addref.overall_coeff).
530                         mul_dyn(ex_to_numeric(overall_coeff))))
531             ->setflag(status_flags::dynallocated  |
532                       status_flags::evaluated );
533     }
534     return this->hold();
535 }
536
537 /*
538 ex mul::eval(int level) const
539 {
540     // simplifications: *(...,x,c1,c2) -> *(...,x,c1*c2) (c1, c2 numeric())
541     //                  *(...,(c1,c2)) -> (...,(c1^c2,1)) (normalize)
542     //                  *(...,x,1) -> +(...,x)
543     //                  *(...,x,0) -> 0
544     //                  *(+(x,y,...),c) -> *(+(*(x,c),*(y,c),...)) (c numeric())
545     //                  *(x) -> x
546     //                  *() -> 1
547
548     debugmsg("mul eval",LOGLEVEL_MEMBER_FUNCTION);
549
550     epvector newseq=seq;
551     epvector::iterator it1,it2;
552     
553     // *(...,x,c1,c2) -> *(...,x,c1*c2) (c1, c2 numeric())
554     it2=newseq.end()-1;
555     it1=it2-1;
556     while ((newseq.size()>=2)&&is_exactly_of_type(*(*it1).rest.bp,numeric)&&
557                                is_exactly_of_type(*(*it2).rest.bp,numeric)) {
558         *it1=expair(ex_to_numeric((*it1).rest).power(ex_to_numeric((*it1).coeff))
559                     .mul(ex_to_numeric((*it2).rest).power(ex_to_numeric((*it2).coeff))),exONE());
560         newseq.pop_back();
561         it2=newseq.end()-1;
562         it1=it2-1;
563     }
564
565     if ((newseq.size()>=1)&&is_exactly_of_type(*(*it2).rest.bp,numeric)) {
566         // *(...,(c1,c2)) -> (...,(c1^c2,1)) (normalize)
567         *it2=expair(ex_to_numeric((*it2).rest).power(ex_to_numeric((*it2).coeff)),exONE());
568         // *(...,x,1) -> *(...,x)
569         if (static_cast<numeric &>(*(*it2).rest.bp).compare(numONE())==0) {
570             newseq.pop_back();
571             it2=newseq.end()-1;
572         }
573     }
574
575     // *(...,x,0) -> 0
576     if ((newseq.size()>=1)&&is_exactly_of_type(*(*it2).rest.bp,numeric)) {
577         if (static_cast<numeric &>(*(*it2).rest.bp).is_zero()==0) {
578             return exZERO();
579         }
580     }
581
582     // *(+(x,y,...),c) -> +(*(x,c),*(y,c),...) (c numeric(), no powers of +())
583     if ((newseq.size()==2)&&is_ex_exactly_of_type(newseq[0].rest,add)&&
584         is_ex_exactly_of_type(newseq[1].rest,numeric)&&
585         (ex_to_numeric(newseq[0].coeff).compare(numONE())==0)) {
586         add const & addref=ex_to_add(newseq[0].rest);
587         numeric const & numref=ex_to_numeric(newseq[1].rest);
588         epvector distrseq;
589         distrseq.reserve(addref.seq.size());
590         for (epvector::const_iterator cit=addref.seq.begin(); cit!=addref.seq.end(); ++cit) {
591             distrseq.push_back(expair((*cit).rest,ex_to_numeric((*cit).coeff).mul(numref)));
592         }
593         return (new add(distrseq,1))->setflag(status_flags::dynallocated  |
594                                               status_flags::evaluated );
595     }
596     
597     if (newseq.size()==0) {
598         // *() -> 1
599         return exONE();
600     } else if (newseq.size()==1) {
601         // *(x) -> x
602         return recombine_pair_to_ex(*newseq.begin());
603     }
604
605     return (new mul(newseq,1))->setflag(status_flags::dynallocated  |
606                                         status_flags::evaluated );
607 }
608 */
609
610 exvector mul::get_indices(void) const
611 {
612     // return union of indices of factors
613     exvector iv;
614     for (epvector::const_iterator cit=seq.begin(); cit!=seq.end(); ++cit) {
615         exvector subiv=(*cit).rest.get_indices();
616         iv.reserve(iv.size()+subiv.size());
617         for (exvector::const_iterator cit2=subiv.begin(); cit2!=subiv.end(); ++cit2) {
618             iv.push_back(*cit2);
619         }
620     }
621     return iv;
622 }
623
624 ex mul::simplify_ncmul(exvector const & v) const
625 {
626     throw(std::logic_error("mul::simplify_ncmul() should never have been called!"));
627 }
628
629 // protected
630
631 int mul::compare_same_type(basic const & other) const
632 {
633     return expairseq::compare_same_type(other);
634 }
635
636 bool mul::is_equal_same_type(basic const & other) const
637 {
638     return expairseq::is_equal_same_type(other);
639 }
640
641 unsigned mul::return_type(void) const
642 {
643     if (seq.size()==0) {
644         // mul without factors: should not happen, but commutes
645         return return_types::commutative;
646     }
647
648     bool all_commutative=1;
649     unsigned rt;
650     epvector::const_iterator cit_noncommutative_element; // point to first found nc element
651
652     for (epvector::const_iterator cit=seq.begin(); cit!=seq.end(); ++cit) {
653         rt=(*cit).rest.return_type();
654         if (rt==return_types::noncommutative_composite) return rt; // one ncc -> mul also ncc
655         if ((rt==return_types::noncommutative)&&(all_commutative)) {
656             // first nc element found, remember position
657             cit_noncommutative_element=cit;
658             all_commutative=0;
659         }
660         if ((rt==return_types::noncommutative)&&(!all_commutative)) {
661                 // another nc element found, compare type_infos
662             if ((*cit_noncommutative_element).rest.return_type_tinfo()!=(*cit).rest.return_type_tinfo()) {
663                 // diffent types -> mul is ncc
664                 return return_types::noncommutative_composite;
665             }
666         }
667     }
668     // all factors checked
669     return all_commutative ? return_types::commutative : return_types::noncommutative;
670 }
671    
672 unsigned mul::return_type_tinfo(void) const
673 {
674     if (seq.size()==0) {
675         // mul without factors: should not happen
676         return tinfo_key;
677     }
678     // return type_info of first noncommutative element
679     for (epvector::const_iterator cit=seq.begin(); cit!=seq.end(); ++cit) {
680         if ((*cit).rest.return_type()==return_types::noncommutative) {
681             return (*cit).rest.return_type_tinfo();
682         }
683     }
684     // no noncommutative element found, should not happen
685     return tinfo_key;
686 }
687
688 ex mul::thisexpairseq(epvector const & v, ex const & oc) const
689 {
690     return (new mul(v,oc))->setflag(status_flags::dynallocated);
691 }
692
693 ex mul::thisexpairseq(epvector * vp, ex const & oc) const
694 {
695     return (new mul(vp,oc))->setflag(status_flags::dynallocated);
696 }
697
698 expair mul::split_ex_to_pair(ex const & e) const
699 {
700     if (is_ex_exactly_of_type(e,power)) {
701         power const & powerref=ex_to_power(e);
702         if (is_ex_exactly_of_type(powerref.exponent,numeric)) {
703             return expair(powerref.basis,powerref.exponent);
704         }
705     }
706     return expair(e,exONE());
707 }
708     
709 expair mul::combine_ex_with_coeff_to_pair(ex const & e,
710                                           ex const & c) const
711 {
712     // to avoid duplication of power simplification rules,
713     // we create a temporary power object
714     // otherwise it would be hard to correctly simplify
715     // expression like (4^(1/3))^(3/2)
716     if (are_ex_trivially_equal(c,exONE())) {
717         return split_ex_to_pair(e);
718     }
719     return split_ex_to_pair(power(e,c));
720 }
721     
722 expair mul::combine_pair_with_coeff_to_pair(expair const & p,
723                                             ex const & c) const
724 {
725     // to avoid duplication of power simplification rules,
726     // we create a temporary power object
727     // otherwise it would be hard to correctly simplify
728     // expression like (4^(1/3))^(3/2)
729     if (are_ex_trivially_equal(c,exONE())) {
730         return p;
731     }
732     return split_ex_to_pair(power(recombine_pair_to_ex(p),c));
733 }
734     
735 ex mul::recombine_pair_to_ex(expair const & p) const
736 {
737     // if (p.coeff.compare(exONE())==0) {
738     // if (are_ex_trivially_equal(p.coeff,exONE())) {
739     if (ex_to_numeric(p.coeff).is_equal(numONE())) {
740         return p.rest;
741     } else {
742         return power(p.rest,p.coeff);
743     }
744 }
745
746 bool mul::expair_needs_further_processing(epp it)
747 {
748     if (is_ex_exactly_of_type((*it).rest,mul) &&
749         ex_to_numeric((*it).coeff).is_integer()) {
750         // combined pair is product with integer power -> expand it
751         *it=split_ex_to_pair(recombine_pair_to_ex(*it));
752         return true;
753     }
754     if (is_ex_exactly_of_type((*it).rest,numeric)) {
755         expair ep=split_ex_to_pair(recombine_pair_to_ex(*it));
756         if (!ep.is_equal(*it)) {
757             // combined pair is a numeric power which can be simplified
758             *it=ep;
759             return true;
760         }
761         if (ex_to_numeric((*it).coeff).is_equal(numONE())) {
762             // combined pair has coeff 1 and must be moved to the end
763             return true;
764         }
765     }
766     return false;
767 }       
768
769 ex mul::default_overall_coeff(void) const
770 {
771     return exONE();
772 }
773
774 void mul::combine_overall_coeff(ex const & c)
775 {
776     ASSERT(is_ex_exactly_of_type(overall_coeff,numeric));
777     ASSERT(is_ex_exactly_of_type(c,numeric));
778     overall_coeff = ex_to_numeric(overall_coeff).mul_dyn(ex_to_numeric(c));
779 }
780
781 void mul::combine_overall_coeff(ex const & c1, ex const & c2)
782 {
783     ASSERT(is_ex_exactly_of_type(overall_coeff,numeric));
784     ASSERT(is_ex_exactly_of_type(c1,numeric));
785     ASSERT(is_ex_exactly_of_type(c2,numeric));
786     overall_coeff = ex_to_numeric(overall_coeff).
787                         mul_dyn(ex_to_numeric(c1).power(ex_to_numeric(c2)));
788 }
789
790 bool mul::can_make_flat(expair const & p) const
791 {
792     ASSERT(is_ex_exactly_of_type(p.coeff,numeric));
793     // this assertion will probably fail somewhere
794     // it would require a more careful make_flat, obeying the power laws
795     // probably should return true only if p.coeff is integer
796     return ex_to_numeric(p.coeff).is_equal(numONE());
797 }
798
799 ex mul::expand(unsigned options) const
800 {
801     exvector sub_expanded_seq;
802     intvector positions_of_adds;
803     intvector number_of_add_operands;
804
805     epvector * expanded_seqp=expandchildren(options);
806
807     epvector const & expanded_seq = expanded_seqp==0 ? seq : *expanded_seqp;
808
809     positions_of_adds.resize(expanded_seq.size());
810     number_of_add_operands.resize(expanded_seq.size());
811
812     int number_of_adds=0;
813     int number_of_expanded_terms=1;
814
815     unsigned current_position=0;
816     epvector::const_iterator last=expanded_seq.end();
817     for (epvector::const_iterator cit=expanded_seq.begin(); cit!=last; ++cit) {
818         if (is_ex_exactly_of_type((*cit).rest,add)&&
819             (ex_to_numeric((*cit).coeff).is_equal(numONE()))) {
820             positions_of_adds[number_of_adds]=current_position;
821             add const & expanded_addref=ex_to_add((*cit).rest);
822             int addref_nops=expanded_addref.nops();
823             number_of_add_operands[number_of_adds]=addref_nops;
824             number_of_expanded_terms *= addref_nops;
825             number_of_adds++;
826         }
827         current_position++;
828     }
829
830     if (number_of_adds==0) {
831         if (expanded_seqp==0) {
832             return this->setflag(status_flags::expanded);
833         }
834         return (new mul(expanded_seqp,overall_coeff))->
835                      setflag(status_flags::dynallocated ||
836                              status_flags::expanded);
837     }
838
839     exvector distrseq;
840     distrseq.reserve(number_of_expanded_terms);
841
842     intvector k;
843     k.resize(number_of_adds);
844     
845     int l;
846     for (l=0; l<number_of_adds; l++) {
847         k[l]=0;
848     }
849
850     while (1) {
851         epvector term;
852         term=expanded_seq;
853         for (l=0; l<number_of_adds; l++) {
854             add const & addref=ex_to_add(expanded_seq[positions_of_adds[l]].rest);
855             ASSERT(term[positions_of_adds[l]].coeff.compare(exONE())==0);
856             term[positions_of_adds[l]]=split_ex_to_pair(addref.op(k[l]));
857         }
858         /*
859         cout << "mul::expand() term begin" << endl;
860         for (epvector::const_iterator cit=term.begin(); cit!=term.end(); ++cit) {
861             cout << "rest" << endl;
862             (*cit).rest.printtree(cout);
863             cout << "coeff" << endl;
864             (*cit).coeff.printtree(cout);
865         }
866         cout << "mul::expand() term end" << endl;
867         */
868         distrseq.push_back((new mul(term,overall_coeff))->
869                                 setflag(status_flags::dynallocated |
870                                         status_flags::expanded));
871
872         // increment k[]
873         l=number_of_adds-1;
874         while ((l>=0)&&((++k[l])>=number_of_add_operands[l])) {
875             k[l]=0;    
876             l--;
877         }
878         if (l<0) break;
879     }
880
881     if (expanded_seqp!=0) {
882         delete expanded_seqp;
883     }
884     /*
885     cout << "mul::expand() distrseq begin" << endl;
886     for (exvector::const_iterator cit=distrseq.begin(); cit!=distrseq.end(); ++cit) {
887         (*cit).printtree(cout);
888     }
889     cout << "mul::expand() distrseq end" << endl;
890     */
891
892     return (new add(distrseq))->setflag(status_flags::dynallocated |
893                                         status_flags::expanded);
894 }
895
896 /*
897 ex mul::expand(unsigned options) const
898 {
899     exvector sub_expanded_seq;
900     intvector positions_of_adds;
901     intvector number_of_add_operands;
902
903     sub_expanded_seq.resize(seq.size());
904     positions_of_adds.resize(seq.size());
905     number_of_add_operands.reserve(seq.size());
906
907     int number_of_adds=0;
908     int number_of_expanded_terms=1;
909     for (unsigned current_position=0; current_position<seq.size(); current_position++) {
910         ex const & expanded_ex=recombine_pair_to_ex(seq[current_position]).expand(options);
911         if (is_ex_exactly_of_type(expanded_ex,add)) {
912             positions_of_adds[number_of_adds]=current_position;
913             add const & expanded_addref=ex_to_add(expanded_ex);
914             number_of_add_operands[number_of_adds]=expanded_addref.seq.size();
915             number_of_expanded_terms *= expanded_addref.seq.size();
916             number_of_adds++;
917         }
918         sub_expanded_seq.push_back(expanded_ex);
919     }
920
921     exvector distrseq;
922     distrseq.reserve(number_of_expanded_terms);
923
924     intvector k;
925     k.resize(number_of_adds);
926     
927     int l;
928     for (l=0; l<number_of_adds; l++) {
929         k[l]=0;
930     }
931
932     while (1) {
933         exvector term;
934         term=sub_expanded_seq;
935         for (l=0; l<number_of_adds; l++) {
936             add const & addref=ex_to_add(sub_expanded_seq[positions_of_adds[l]]);
937             term[positions_of_adds[l]]=addref.recombine_pair_to_ex(addref.seq[k[l]]);
938         }
939         distrseq.push_back((new mul(term))->setflag(status_flags::dynallocated |
940                                                     status_flags::expanded));
941
942         // increment k[]
943         l=number_of_adds-1;
944         while ((l>=0)&&((++k[l])>=number_of_add_operands[l])) {
945             k[l]=0;    
946             l--;
947         }
948         if (l<0) break;
949     }
950
951     return (new add(distrseq))->setflag(status_flags::dynallocated |
952                                         status_flags::expanded);
953 }
954 */
955
956 //////////
957 // new virtual functions which can be overridden by derived classes
958 //////////
959
960 // none
961
962 //////////
963 // non-virtual functions in this class
964 //////////
965
966 epvector * mul::expandchildren(unsigned options) const
967 {
968     epvector::const_iterator last=seq.end();
969     epvector::const_iterator cit=seq.begin();
970     while (cit!=last) {
971         ex const & factor=recombine_pair_to_ex(*cit);
972         ex const & expanded_factor=factor.expand(options);
973         if (!are_ex_trivially_equal(factor,expanded_factor)) {
974
975             // something changed, copy seq, eval and return it
976             epvector *s=new epvector;
977             s->reserve(seq.size());
978
979             // copy parts of seq which are known not to have changed
980             epvector::const_iterator cit2=seq.begin();
981             while (cit2!=cit) {
982                 s->push_back(*cit2);
983                 ++cit2;
984             }
985             // copy first changed element
986             s->push_back(split_ex_to_pair(expanded_factor));
987             ++cit2;
988             // copy rest
989             while (cit2!=last) {
990                 s->push_back(split_ex_to_pair(recombine_pair_to_ex(*cit2).expand(options)));
991                 ++cit2;
992             }
993             return s;
994         }
995         ++cit;
996     }
997     
998     return 0; // nothing has changed
999 }
1000    
1001 //////////
1002 // static member variables
1003 //////////
1004
1005 // protected
1006
1007 unsigned mul::precedence=50;
1008
1009
1010 //////////
1011 // global constants
1012 //////////
1013
1014 const mul some_mul;
1015 type_info const & typeid_mul=typeid(some_mul);
1016
1017 } // namespace GiNaC