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