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