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