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