- switched to automake build environment
[ginac.git] / ginac / ncmul.cpp
1 /** @file ncmul.cpp
2  *
3  *  Implementation of GiNaC's non-commutative products of expressions. */
4
5 #include <algorithm>
6 #include <iostream>
7 #include <stdexcept>
8
9 #include "ginac.h"
10
11 //////////
12 // default constructor, destructor, copy constructor assignment operator and helpers
13 //////////
14
15 // public
16
17 ncmul::ncmul()
18 {
19     debugmsg("ncmul default constructor",LOGLEVEL_CONSTRUCT);
20     tinfo_key = TINFO_NCMUL;
21 }
22
23 ncmul::~ncmul()
24 {
25     debugmsg("ncmul destructor",LOGLEVEL_DESTRUCT);
26     destroy(0);
27 }
28
29 ncmul::ncmul(ncmul const & other)
30 {
31     debugmsg("ncmul copy constructor",LOGLEVEL_CONSTRUCT);
32     copy(other);
33 }
34
35 ncmul const & ncmul::operator=(ncmul const & other)
36 {
37     debugmsg("ncmul operator=",LOGLEVEL_ASSIGNMENT);
38     if (this != &other) {
39         destroy(1);
40         copy(other);
41     }
42     return *this;
43 }
44
45 // protected
46
47 void ncmul::copy(ncmul const & other)
48 {
49     exprseq::copy(other);
50 }
51
52 void ncmul::destroy(bool call_parent)
53 {
54     if (call_parent) exprseq::destroy(call_parent);
55 }
56
57 //////////
58 // other constructors
59 //////////
60
61 // public
62
63 ncmul::ncmul(ex const & lh, ex const & rh) :
64     exprseq(lh,rh)
65 {
66     debugmsg("ncmul constructor from ex,ex",LOGLEVEL_CONSTRUCT);
67     tinfo_key = TINFO_NCMUL;
68 }
69
70 ncmul::ncmul(ex const & f1, ex const & f2, ex const & f3) :
71     exprseq(f1,f2,f3)
72 {
73     debugmsg("ncmul constructor from 3 ex",LOGLEVEL_CONSTRUCT);
74     tinfo_key = TINFO_NCMUL;
75 }
76
77 ncmul::ncmul(ex const & f1, ex const & f2, ex const & f3,
78       ex const & f4) : exprseq(f1,f2,f3,f4)
79 {
80     debugmsg("ncmul constructor from 4 ex",LOGLEVEL_CONSTRUCT);
81     tinfo_key = TINFO_NCMUL;
82 }
83
84 ncmul::ncmul(ex const & f1, ex const & f2, ex const & f3,
85       ex const & f4, ex const & f5) : exprseq(f1,f2,f3,f4,f5)
86 {
87     debugmsg("ncmul constructor from 5 ex",LOGLEVEL_CONSTRUCT);
88     tinfo_key = TINFO_NCMUL;
89 }
90
91 ncmul::ncmul(ex const & f1, ex const & f2, ex const & f3,
92       ex const & f4, ex const & f5, ex const & f6) :
93     exprseq(f1,f2,f3,f4,f5,f6)
94 {
95     debugmsg("ncmul constructor from 6 ex",LOGLEVEL_CONSTRUCT);
96     tinfo_key = TINFO_NCMUL;
97 }
98
99 ncmul::ncmul(exvector const & v, bool discardable) : exprseq(v,discardable)
100 {
101     debugmsg("ncmul constructor from exvector,bool",LOGLEVEL_CONSTRUCT);
102     tinfo_key = TINFO_NCMUL;
103 }
104
105 ncmul::ncmul(exvector * vp) : exprseq(vp)
106 {
107     debugmsg("ncmul constructor from exvector *",LOGLEVEL_CONSTRUCT);
108     tinfo_key = TINFO_NCMUL;
109 }
110     
111 //////////
112 // functions overriding virtual functions from bases classes
113 //////////
114
115 // public
116
117 basic * ncmul::duplicate() const
118 {
119     debugmsg("ncmul duplicate",LOGLEVEL_ASSIGNMENT);
120     return new ncmul(*this);
121 }
122
123 bool ncmul::info(unsigned inf) const
124 {
125     throw(std::logic_error("which flags have to be implemented in ncmul::info()?"));
126 }
127
128 typedef vector<int> intvector;
129
130 ex ncmul::expand(unsigned options) const
131 {
132     exvector sub_expanded_seq;
133     intvector positions_of_adds;
134     intvector number_of_add_operands;
135
136     exvector expanded_seq=expandchildren(options);
137
138     positions_of_adds.resize(expanded_seq.size());
139     number_of_add_operands.resize(expanded_seq.size());
140
141     int number_of_adds=0;
142     int number_of_expanded_terms=1;
143
144     unsigned current_position=0;
145     exvector::const_iterator last=expanded_seq.end();
146     for (exvector::const_iterator cit=expanded_seq.begin(); cit!=last; ++cit) {
147         if (is_ex_exactly_of_type((*cit),add)) {
148             positions_of_adds[number_of_adds]=current_position;
149             add const & expanded_addref=ex_to_add(*cit);
150             number_of_add_operands[number_of_adds]=expanded_addref.seq.size();
151             number_of_expanded_terms *= expanded_addref.seq.size();
152             number_of_adds++;
153         }
154         current_position++;
155     }
156
157     if (number_of_adds==0) {
158         return (new ncmul(expanded_seq,1))->setflag(status_flags::dynallocated ||
159                                                     status_flags::expanded);
160     }
161
162     exvector distrseq;
163     distrseq.reserve(number_of_expanded_terms);
164
165     intvector k;
166     k.resize(number_of_adds);
167     
168     int l;
169     for (l=0; l<number_of_adds; l++) {
170         k[l]=0;
171     }
172
173     while (1) {
174         exvector term;
175         term=expanded_seq;
176         for (l=0; l<number_of_adds; l++) {
177             ASSERT(is_ex_exactly_of_type(expanded_seq[positions_of_adds[l]],add));
178             add const & addref=ex_to_add(expanded_seq[positions_of_adds[l]]);
179             term[positions_of_adds[l]]=addref.recombine_pair_to_ex(addref.seq[k[l]]);
180         }
181         distrseq.push_back((new ncmul(term,1))->setflag(status_flags::dynallocated |
182                                                         status_flags::expanded));
183
184         // increment k[]
185         l=number_of_adds-1;
186         while ((l>=0)&&((++k[l])>=number_of_add_operands[l])) {
187             k[l]=0;    
188             l--;
189         }
190         if (l<0) break;
191     }
192
193     return (new add(distrseq))->setflag(status_flags::dynallocated |
194                                         status_flags::expanded);
195 }
196
197 int ncmul::degree(symbol const & s) const
198 {
199     int deg_sum=0;
200     for (exvector::const_iterator cit=seq.begin(); cit!=seq.end(); ++cit) {
201         deg_sum+=(*cit).degree(s);
202     }
203     return deg_sum;
204 }
205
206 int ncmul::ldegree(symbol const & s) const
207 {
208     int deg_sum=0;
209     for (exvector::const_iterator cit=seq.begin(); cit!=seq.end(); ++cit) {
210         deg_sum+=(*cit).ldegree(s);
211     }
212     return deg_sum;
213 }
214
215 ex ncmul::coeff(symbol const & s, int const n) const
216 {
217     exvector coeffseq;
218     coeffseq.reserve(seq.size());
219
220     if (n==0) {
221         // product of individual coeffs
222         // if a non-zero power of s is found, the resulting product will be 0
223         exvector::const_iterator it=seq.begin();
224         while (it!=seq.end()) {
225             coeffseq.push_back((*it).coeff(s,n));
226             ++it;
227         }
228         return (new ncmul(coeffseq,1))->setflag(status_flags::dynallocated);
229     }
230          
231     exvector::const_iterator it=seq.begin();
232     bool coeff_found=0;
233     while (it!=seq.end()) {
234         ex c=(*it).coeff(s,n);
235         if (!c.is_zero()) {
236             coeffseq.push_back(c);
237             coeff_found=1;
238         } else {
239             coeffseq.push_back(*it);
240         }
241         ++it;
242     }
243
244     if (coeff_found) return (new ncmul(coeffseq,1))->setflag(status_flags::dynallocated);
245     
246     return exZERO();
247 }
248
249 unsigned ncmul::count_factors(ex const & e) const
250 {
251     if ((is_ex_exactly_of_type(e,mul)&&(e.return_type()!=return_types::commutative))||
252         (is_ex_exactly_of_type(e,ncmul))) {
253         unsigned factors=0;
254         for (int i=0; i<e.nops(); i++) {
255             factors += count_factors(e.op(i));
256         }
257         return factors;
258     }
259     return 1;
260 }
261         
262 void ncmul::append_factors(exvector & v, ex const & e) const
263 {
264     if ((is_ex_exactly_of_type(e,mul)&&(e.return_type()!=return_types::commutative))||
265         (is_ex_exactly_of_type(e,ncmul))) {
266         for (int i=0; i<e.nops(); i++) {
267             append_factors(v,e.op(i));
268         }
269         return;
270     }
271     v.push_back(e);
272 }
273
274 typedef vector<unsigned> unsignedvector;
275 typedef vector<exvector> exvectorvector;
276
277 ex ncmul::eval(int level) const
278 {
279     // simplifications: ncmul(...,*(x1,x2),...,ncmul(x3,x4),...) ->
280     //                      ncmul(...,x1,x2,...,x3,x4,...) (associativity)
281     //                  ncmul(x) -> x
282     //                  ncmul() -> 1
283     //                  ncmul(...,c1,...,c2,...) ->
284     //                      *(c1,c2,ncmul(...)) (pull out commutative elements)
285     //                  ncmul(x1,y1,x2,y2) -> *(ncmul(x1,x2),ncmul(y1,y2))
286     //                      (collect elements of same type)
287     //                  ncmul(x1,x2,x3,...) -> x::eval_ncmul(x1,x2,x3,...)
288     // the following rule would be nice, but produces a recursion,
289     // which must be trapped by introducing a flag that the sub-ncmuls()
290     // are already evaluated (maybe later...)
291     //                  ncmul(x1,x2,...,X,y1,y2,...) ->
292     //                      ncmul(ncmul(x1,x2,...),X,ncmul(y1,y2,...)
293     //                      (X noncommutative_composite)
294
295     if ((level==1)&&(flags & status_flags::evaluated)) {
296         return *this;
297     }
298
299     exvector evaledseq=evalchildren(level);
300
301     // ncmul(...,*(x1,x2),...,ncmul(x3,x4),...) ->
302     //     ncmul(...,x1,x2,...,x3,x4,...) (associativity)
303     unsigned factors=0;
304     for (exvector::const_iterator cit=evaledseq.begin(); cit!=evaledseq.end(); ++cit) {
305         factors += count_factors(*cit);
306     }
307
308     exvector assocseq;
309     assocseq.reserve(factors);
310     for (exvector::const_iterator cit=evaledseq.begin(); cit!=evaledseq.end(); ++cit) {
311         append_factors(assocseq,*cit);
312     }
313
314     // ncmul(x) -> x
315     if (assocseq.size()==1) return *(seq.begin());
316
317     // ncmul() -> 1
318     if (assocseq.size()==0) return exONE();
319
320     // determine return types
321     unsignedvector rettypes;
322     rettypes.reserve(assocseq.size());
323     unsigned i=0;
324     unsigned count_commutative=0;
325     unsigned count_noncommutative=0;
326     unsigned count_noncommutative_composite=0;
327     for (exvector::const_iterator cit=assocseq.begin(); cit!=assocseq.end(); ++cit) {
328         switch (rettypes[i]=(*cit).return_type()) {
329         case return_types::commutative:
330             count_commutative++;
331             break;
332         case return_types::noncommutative:
333             count_noncommutative++;
334             break;
335         case return_types::noncommutative_composite:
336             count_noncommutative_composite++;
337             break;
338         default:
339             throw(std::logic_error("ncmul::eval(): invalid return type"));
340         }
341         ++i;
342     }
343     ASSERT(count_commutative+count_noncommutative+count_noncommutative_composite==assocseq.size());
344
345     // ncmul(...,c1,...,c2,...) ->
346     //     *(c1,c2,ncmul(...)) (pull out commutative elements)
347     if (count_commutative!=0) {
348         exvector commutativeseq;
349         commutativeseq.reserve(count_commutative+1);
350         exvector noncommutativeseq;
351         noncommutativeseq.reserve(assocseq.size()-count_commutative);
352         for (i=0; i<assocseq.size(); ++i) {
353             if (rettypes[i]==return_types::commutative) {
354                 commutativeseq.push_back(assocseq[i]);
355             } else {
356                 noncommutativeseq.push_back(assocseq[i]);
357             }
358         }
359         commutativeseq.push_back((new ncmul(noncommutativeseq,1))->
360                                   setflag(status_flags::dynallocated));
361         return (new mul(commutativeseq))->setflag(status_flags::dynallocated);
362     }
363         
364     // ncmul(x1,y1,x2,y2) -> *(ncmul(x1,x2),ncmul(y1,y2))
365     //     (collect elements of same type)
366
367     if (count_noncommutative_composite==0) {
368         // there are neither commutative nor noncommutative_composite
369         // elements in assocseq
370         ASSERT(count_commutative==0);
371
372         exvectorvector evv;
373         unsignedvector rttinfos;
374         evv.reserve(assocseq.size());
375         rttinfos.reserve(assocseq.size());
376
377         for (exvector::const_iterator cit=assocseq.begin(); cit!=assocseq.end(); ++cit) {
378             unsigned ti=(*cit).return_type_tinfo();
379             // search type in vector of known types
380             for (i=0; i<rttinfos.size(); ++i) {
381                 if (ti==rttinfos[i]) {
382                     evv[i].push_back(*cit);
383                     break;
384                 }
385             }
386             if (i>=rttinfos.size()) {
387                 // new type
388                 rttinfos.push_back(ti);
389                 evv.push_back(exvector());
390                 (*(evv.end()-1)).reserve(assocseq.size());
391                 (*(evv.end()-1)).push_back(*cit);
392             }
393         }
394
395 #ifdef DOASSERT
396         ASSERT(evv.size()==rttinfos.size());
397         ASSERT(evv.size()>0);
398         unsigned s=0;
399         for (i=0; i<evv.size(); ++i) {
400             s += evv[i].size();
401         }
402         ASSERT(s==assocseq.size());
403 #endif // def DOASSERT
404         
405         // if all elements are of same type, simplify the string
406         if (evv.size()==1) {
407             return evv[0][0].simplify_ncmul(evv[0]);
408         }
409         
410         exvector splitseq;
411         splitseq.reserve(evv.size());
412         for (i=0; i<evv.size(); ++i) {
413             splitseq.push_back((new ncmul(evv[i]))->
414                                setflag(status_flags::dynallocated));
415         }
416
417         return (new mul(splitseq))->setflag(status_flags::dynallocated);
418     }
419     
420     return (new ncmul(assocseq))->setflag(status_flags::dynallocated |
421                                           status_flags::evaluated);
422 }
423
424 exvector ncmul::get_indices(void) const
425 {
426     // return union of indices of factors
427     exvector iv;
428     for (exvector::const_iterator cit=seq.begin(); cit!=seq.end(); ++cit) {
429         exvector subiv=(*cit).get_indices();
430         iv.reserve(iv.size()+subiv.size());
431         for (exvector::const_iterator cit2=subiv.begin(); cit2!=subiv.end(); ++cit2) {
432             iv.push_back(*cit2);
433         }
434     }
435     return iv;
436 }
437
438 ex ncmul::subs(lst const & ls, lst const & lr) const
439 {
440     return ncmul(subschildren(ls, lr));
441 }
442
443 ex ncmul::thisexprseq(exvector const & v) const
444 {
445     return (new ncmul(v))->setflag(status_flags::dynallocated);
446 }
447
448 ex ncmul::thisexprseq(exvector * vp) const
449 {
450     return (new ncmul(vp))->setflag(status_flags::dynallocated);
451 }
452
453 // protected
454
455 int ncmul::compare_same_type(basic const & other) const
456 {
457     return exprseq::compare_same_type(other);
458 }
459
460 unsigned ncmul::return_type(void) const
461 {
462     if (seq.size()==0) {
463         // ncmul without factors: should not happen, but commutes
464         return return_types::commutative;
465     }
466
467     bool all_commutative=1;
468     unsigned rt;
469     exvector::const_iterator cit_noncommutative_element; // point to first found nc element
470
471     for (exvector::const_iterator cit=seq.begin(); cit!=seq.end(); ++cit) {
472         rt=(*cit).return_type();
473         if (rt==return_types::noncommutative_composite) return rt; // one ncc -> mul also ncc
474         if ((rt==return_types::noncommutative)&&(all_commutative)) {
475             // first nc element found, remember position
476             cit_noncommutative_element=cit;
477             all_commutative=0;
478         }
479         if ((rt==return_types::noncommutative)&&(!all_commutative)) {
480             // another nc element found, compare type_infos
481             if ((*cit_noncommutative_element).return_type_tinfo()!=(*cit).return_type_tinfo()) {
482                 // diffent types -> mul is ncc
483                 return return_types::noncommutative_composite;
484             }
485         }
486     }
487     // all factors checked
488     ASSERT(!all_commutative); // not all factors should commute, because this is a ncmul();
489     return all_commutative ? return_types::commutative : return_types::noncommutative;
490 }
491    
492 unsigned ncmul::return_type_tinfo(void) const
493 {
494     if (seq.size()==0) {
495         // mul without factors: should not happen
496         return tinfo_key;
497     }
498     // return type_info of first noncommutative element
499     for (exvector::const_iterator cit=seq.begin(); cit!=seq.end(); ++cit) {
500         if ((*cit).return_type()==return_types::noncommutative) {
501             return (*cit).return_type_tinfo();
502         }
503     }
504     // no noncommutative element found, should not happen
505     return tinfo_key;
506 }
507
508 //////////
509 // new virtual functions which can be overridden by derived classes
510 //////////
511
512 // none
513
514 //////////
515 // non-virtual functions in this class
516 //////////
517
518 exvector ncmul::expandchildren(unsigned options) const
519 {
520     exvector s;
521     s.reserve(seq.size());
522
523     for (exvector::const_iterator it=seq.begin(); it!=seq.end(); ++it) {
524         s.push_back((*it).expand(options));
525     }
526     return s;
527 }
528
529 exvector const & ncmul::get_factors(void) const
530 {
531     return seq;
532 }
533
534 //////////
535 // static member variables
536 //////////
537
538 // protected
539
540 unsigned ncmul::precedence=50;
541
542
543 //////////
544 // global constants
545 //////////
546
547 const ncmul some_ncmul;
548 type_info const & typeid_ncmul=typeid(some_ncmul);
549
550 //////////
551 // friend functions
552 //////////
553
554 ex nonsimplified_ncmul(exvector const & v)
555 {
556     return (new ncmul(v))->setflag(status_flags::dynallocated);
557 }
558
559 ex simplified_ncmul(exvector const & v)
560 {
561     if (v.size()==0) {
562         return exONE();
563     } else if (v.size()==1) {
564         return v[0];
565     }
566     return (new ncmul(v))->setflag(status_flags::dynallocated |
567                                    status_flags::evaluated);
568 }
569
570