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