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