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