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