- ASSERT macro renamed to GINAC_ASSERT
[ginac.git] / ginac / add.cpp
1 /** @file add.cpp
2  *
3  *  Implementation of GiNaC's sums 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 <iostream>
24 #include <stdexcept>
25
26 #include "add.h"
27 #include "mul.h"
28 #include "debugmsg.h"
29
30 namespace GiNaC {
31
32 //////////
33 // default constructor, destructor, copy constructor assignment operator and helpers
34 //////////
35
36 // public
37
38 add::add()
39 {
40     debugmsg("add default constructor",LOGLEVEL_CONSTRUCT);
41     tinfo_key = TINFO_add;
42 }
43
44 add::~add()
45 {
46     debugmsg("add destructor",LOGLEVEL_DESTRUCT);
47     destroy(0);
48 }
49
50 add::add(add const & other)
51 {
52     debugmsg("add copy constructor",LOGLEVEL_CONSTRUCT);
53     copy(other);
54 }
55
56 add const & add::operator=(add const & other)
57 {
58     debugmsg("add operator=",LOGLEVEL_ASSIGNMENT);
59     if (this != &other) {
60         destroy(1);
61         copy(other);
62     }
63     return *this;
64 }
65
66 // protected
67
68 void add::copy(add const & other)
69 {
70     expairseq::copy(other);
71 }
72
73 void add::destroy(bool call_parent)
74 {
75     if (call_parent) expairseq::destroy(call_parent);
76 }
77
78 //////////
79 // other constructors
80 //////////
81
82 // public
83
84 add::add(ex const & lh, ex const & rh)
85 {
86     debugmsg("add constructor from ex,ex",LOGLEVEL_CONSTRUCT);
87     tinfo_key = TINFO_add;
88     overall_coeff=exZERO();
89     construct_from_2_ex(lh,rh);
90     GINAC_ASSERT(is_canonical());
91 }
92
93 add::add(exvector const & v)
94 {
95     debugmsg("add constructor from exvector",LOGLEVEL_CONSTRUCT);
96     tinfo_key = TINFO_add;
97     overall_coeff=exZERO();
98     construct_from_exvector(v);
99     GINAC_ASSERT(is_canonical());
100 }
101
102 /*
103 add::add(epvector const & v, bool do_not_canonicalize)
104 {
105     debugmsg("add constructor from epvector,bool",LOGLEVEL_CONSTRUCT);
106     tinfo_key = TINFO_add;
107     if (do_not_canonicalize) {
108         seq=v;
109 #ifdef EXPAIRSEQ_USE_HASHTAB
110         combine_same_terms(); // to build hashtab
111 #endif // def EXPAIRSEQ_USE_HASHTAB
112     } else {
113         construct_from_epvector(v);
114     }
115     GINAC_ASSERT(is_canonical());
116 }
117 */
118
119 add::add(epvector const & v)
120 {
121     debugmsg("add constructor from epvector",LOGLEVEL_CONSTRUCT);
122     tinfo_key = TINFO_add;
123     overall_coeff=exZERO();
124     construct_from_epvector(v);
125     GINAC_ASSERT(is_canonical());
126 }
127
128 add::add(epvector const & v, ex const & oc)
129 {
130     debugmsg("add constructor from epvector,ex",LOGLEVEL_CONSTRUCT);
131     tinfo_key = TINFO_add;
132     overall_coeff=oc;
133     construct_from_epvector(v);
134     GINAC_ASSERT(is_canonical());
135 }
136
137 add::add(epvector * vp, ex const & oc)
138 {
139     debugmsg("add constructor from epvector *,ex",LOGLEVEL_CONSTRUCT);
140     tinfo_key = TINFO_add;
141     GINAC_ASSERT(vp!=0);
142     overall_coeff=oc;
143     construct_from_epvector(*vp);
144     delete vp;
145     GINAC_ASSERT(is_canonical());
146 }
147
148 //////////
149 // functions overriding virtual functions from bases classes
150 //////////
151
152 // public
153
154 basic * add::duplicate() const
155 {
156     debugmsg("add duplicate",LOGLEVEL_DUPLICATE);
157     return new add(*this);
158 }
159
160 bool add::info(unsigned inf) const
161 {
162     // TODO: optimize
163     if (inf==info_flags::polynomial || inf==info_flags::integer_polynomial || inf==info_flags::rational_polynomial || inf==info_flags::rational_function) {
164         for (epvector::const_iterator it=seq.begin(); it!=seq.end(); ++it) {
165             if (!(recombine_pair_to_ex(*it).info(inf)))
166                 return false;
167         }
168         return true;
169     } else {
170         return expairseq::info(inf);
171     }
172 }
173
174 int add::degree(symbol const & s) const
175 {
176     int deg=INT_MIN;
177     if (!overall_coeff.is_equal(exZERO())) {
178         deg=0;
179     }
180     int cur_deg;
181     for (epvector::const_iterator cit=seq.begin(); cit!=seq.end(); ++cit) {
182         cur_deg=(*cit).rest.degree(s);
183         if (cur_deg>deg) deg=cur_deg;
184     }
185     return deg;
186 }
187
188 int add::ldegree(symbol const & s) const
189 {
190     int deg=INT_MAX;
191     if (!overall_coeff.is_equal(exZERO())) {
192         deg=0;
193     }
194     int cur_deg;
195     for (epvector::const_iterator cit=seq.begin(); cit!=seq.end(); ++cit) {
196         cur_deg=(*cit).rest.ldegree(s);
197         if (cur_deg<deg) deg=cur_deg;
198     }
199     return deg;
200 }
201
202 ex add::coeff(symbol const & s, int const n) const
203 {
204     epvector coeffseq;
205     coeffseq.reserve(seq.size());
206
207     epvector::const_iterator it=seq.begin();
208     while (it!=seq.end()) {
209         coeffseq.push_back(combine_ex_with_coeff_to_pair((*it).rest.coeff(s,n),
210                                                          (*it).coeff));
211         ++it;
212     }
213     if (n==0) {
214         return (new add(coeffseq,overall_coeff))->setflag(status_flags::dynallocated);
215     }
216     return (new add(coeffseq))->setflag(status_flags::dynallocated);
217 }
218
219 /*
220 ex add::eval(int level) const
221 {
222     // simplifications: +(...,x,c1,c2) -> +(...,x,c1+c2) (c1, c2 numeric())
223     //                  +(...,(c1,c2)) -> (...,(c1*c2,1)) (normalize)
224     //                  +(...,x,0) -> +(...,x)
225     //                  +(x) -> x
226     //                  +() -> 0
227
228     debugmsg("add eval",LOGLEVEL_MEMBER_FUNCTION);
229
230     epvector newseq=seq;
231     epvector::iterator it1,it2;
232     
233     // +(...,x,c1,c2) -> +(...,x,c1+c2) (c1, c2 numeric())
234     it2=newseq.end()-1;
235     it1=it2-1;
236     while ((newseq.size()>=2)&&is_exactly_of_type(*(*it1).rest.bp,numeric)&&
237                                is_exactly_of_type(*(*it2).rest.bp,numeric)) {
238         *it1=expair(ex_to_numeric((*it1).rest).mul(ex_to_numeric((*it1).coeff))
239                     .add(ex_to_numeric((*it2).rest).mul(ex_to_numeric((*it2).coeff))),exONE());
240         newseq.pop_back();
241         it2=newseq.end()-1;
242         it1=it2-1;
243     }
244
245     if ((newseq.size()>=1)&&is_exactly_of_type(*(*it2).rest.bp,numeric)) {
246         // +(...,(c1,c2)) -> (...,(c1*c2,1)) (normalize)
247         *it2=expair(ex_to_numeric((*it2).rest).mul(ex_to_numeric((*it2).coeff)),exONE());
248         // +(...,x,0) -> +(...,x)
249         if (ex_to_numeric((*it2).rest).compare(0)==0) {
250             newseq.pop_back();
251         }
252     }
253
254     if (newseq.size()==0) {
255         // +() -> 0
256         return exZERO();
257     } else if (newseq.size()==1) {
258         // +(x) -> x
259         return recombine_pair_to_ex(*(newseq.begin()));
260     }
261
262     return (new add(newseq,1))->setflag(status_flags::dynallocated  |
263                                         status_flags::evaluated );
264 }
265 */
266
267 /*
268 ex add::eval(int level) const
269 {
270     // simplifications: +(...,x,c1,c2) -> +(...,x,c1+c2) (c1, c2 numeric())
271     //                  +(...,(c1,c2)) -> (...,(c1*c2,1)) (normalize)
272     //                  +(...,x,0) -> +(...,x)
273     //                  +(x) -> x
274     //                  +() -> 0
275
276     debugmsg("add eval",LOGLEVEL_MEMBER_FUNCTION);
277
278     if ((level==1)&&(flags & status_flags::evaluated)) {
279 #ifdef DO_GINAC_ASSERT
280         for (epvector::const_iterator cit=seq.begin(); cit!=seq.end(); ++cit) {
281             GINAC_ASSERT(!is_ex_exactly_of_type((*cit).rest,add));
282             GINAC_ASSERT(!(is_ex_exactly_of_type((*cit).rest,numeric)&&
283                      (ex_to_numeric((*cit).coeff).compare(numONE())!=0)));
284         }
285 #endif // def DO_GINAC_ASSERT
286         return *this;
287     }
288
289     epvector newseq;
290     epvector::iterator it1,it2;
291     bool seq_copied=false;
292
293     epvector * evaled_seqp=evalchildren(level);
294     if (evaled_seqp!=0) {
295         // do more evaluation later
296         return (new add(evaled_seqp))->setflag(status_flags::dynallocated);
297     }
298
299 #ifdef DO_GINAC_ASSERT
300     for (epvector::const_iterator cit=seq.begin(); cit!=seq.end(); ++cit) {
301         GINAC_ASSERT(!is_ex_exactly_of_type((*cit).rest,add));
302         GINAC_ASSERT(!(is_ex_exactly_of_type((*cit).rest,numeric)&&
303                  (ex_to_numeric((*cit).coeff).compare(numONE())!=0)));
304     }
305 #endif // def DO_GINAC_ASSERT
306
307     if (flags & status_flags::evaluated) {
308         return *this;
309     }
310     
311     expair const & last_expair=*(seq.end()-1);
312     expair const & next_to_last_expair=*(seq.end()-2);
313     int seq_size = seq.size();
314
315     // +(...,x,c1,c2) -> +(...,x,c1+c2) (c1, c2 numeric())
316     if ((!seq_copied)&&(seq_size>=2)&&
317         is_ex_exactly_of_type(last_expair.rest,numeric)&&
318         is_ex_exactly_of_type(next_to_last_expair.rest,numeric)) {
319         newseq=seq;
320         seq_copied=true;
321         it2=newseq.end()-1;
322         it1=it2-1;
323     }
324     while (seq_copied&&(newseq.size()>=2)&&
325            is_ex_exactly_of_type((*it1).rest,numeric)&&
326            is_ex_exactly_of_type((*it2).rest,numeric)) {
327         *it1=expair(ex_to_numeric((*it1).rest).mul(ex_to_numeric((*it1).coeff))
328                     .add_dyn(ex_to_numeric((*it2).rest).mul(ex_to_numeric((*it2).coeff))),exONE());
329         newseq.pop_back();
330         it2=newseq.end()-1;
331         it1=it2-1;
332     }
333
334     // +(...,(c1,c2)) -> (...,(c1*c2,1)) (normalize)
335     if ((!seq_copied)&&(seq_size>=1)&&
336         (is_ex_exactly_of_type(last_expair.rest,numeric))&&
337         (ex_to_numeric(last_expair.coeff).compare(numONE())!=0)) {
338         newseq=seq;
339         seq_copied=true;
340         it2=newseq.end()-1;
341     }
342     if (seq_copied&&(newseq.size()>=1)&&
343         (is_ex_exactly_of_type((*it2).rest,numeric))&&
344         (ex_to_numeric((*it2).coeff).compare(numONE())!=0)) {
345         *it2=expair(ex_to_numeric((*it2).rest).mul_dyn(ex_to_numeric((*it2).coeff)),exONE());
346     }
347         
348     // +(...,x,0) -> +(...,x)
349     if ((!seq_copied)&&(seq_size>=1)&&
350         (is_ex_exactly_of_type(last_expair.rest,numeric))&&
351         (ex_to_numeric(last_expair.rest).is_zero())) {
352         newseq=seq;
353         seq_copied=true;
354         it2=newseq.end()-1;
355     }
356     if (seq_copied&&(newseq.size()>=1)&&
357         (is_ex_exactly_of_type((*it2).rest,numeric))&&
358         (ex_to_numeric((*it2).rest).is_zero())) {
359         newseq.pop_back();
360     }
361
362     // +() -> 0
363     if ((!seq_copied)&&(seq_size==0)) {
364         return exZERO();
365     } else if (seq_copied&&(newseq.size()==0)) {
366         return exZERO();
367     }
368
369     // +(x) -> x
370     if ((!seq_copied)&&(seq_size==1)) {
371         return recombine_pair_to_ex(*(seq.begin()));
372     } else if (seq_copied&&(newseq.size()==1)) {
373         return recombine_pair_to_ex(*(newseq.begin()));
374     }
375
376     if (!seq_copied) return this->hold();
377
378     return (new add(newseq,1))->setflag(status_flags::dynallocated  |
379                                         status_flags::evaluated );
380 }
381 */
382
383 ex add::eval(int level) const
384 {
385     // simplifications: +(;c) -> c
386     //                  +(x;1) -> x
387
388     debugmsg("add eval",LOGLEVEL_MEMBER_FUNCTION);
389
390     epvector * evaled_seqp=evalchildren(level);
391     if (evaled_seqp!=0) {
392         // do more evaluation later
393         return (new add(evaled_seqp,overall_coeff))->
394                    setflag(status_flags::dynallocated);
395     }
396
397 #ifdef DO_GINAC_ASSERT
398     for (epvector::const_iterator cit=seq.begin(); cit!=seq.end(); ++cit) {
399         GINAC_ASSERT(!is_ex_exactly_of_type((*cit).rest,add));
400         if (is_ex_exactly_of_type((*cit).rest,numeric)) {
401             dbgprint();
402         }
403         GINAC_ASSERT(!is_ex_exactly_of_type((*cit).rest,numeric));
404     }
405 #endif // def DO_GINAC_ASSERT
406
407     if (flags & status_flags::evaluated) {
408         GINAC_ASSERT(seq.size()>0);
409         GINAC_ASSERT((seq.size()>1)||!overall_coeff.is_equal(exZERO()));
410         return *this;
411     }
412
413     int seq_size=seq.size();
414     if (seq_size==0) {
415         // +(;c) -> c
416         return overall_coeff;
417     } else if ((seq_size==1)&&overall_coeff.is_equal(exZERO())) {
418         // +(x;0) -> x
419         return recombine_pair_to_ex(*(seq.begin()));
420     }
421     return this->hold();
422 }
423
424 exvector add::get_indices(void) const
425 {
426     // FIXME: all terms in the sum should have the same indices (compatible
427     // tensors) however this is not checked, since there is no function yet
428     // which compares indices (idxvector can be unsorted)
429     if (seq.size()==0) {
430         return exvector();
431     }
432     return (seq.begin())->rest.get_indices();
433 }    
434
435 ex add::simplify_ncmul(exvector const & v) const
436 {
437     if (seq.size()==0) {
438         return expairseq::simplify_ncmul(v);
439     }
440     return (*seq.begin()).rest.simplify_ncmul(v);
441 }    
442
443 // protected
444
445 int add::compare_same_type(basic const & other) const
446 {
447     return expairseq::compare_same_type(other);
448 }
449
450 bool add::is_equal_same_type(basic const & other) const
451 {
452     return expairseq::is_equal_same_type(other);
453 }
454
455 unsigned add::return_type(void) const
456 {
457     if (seq.size()==0) {
458         return return_types::commutative;
459     }
460     return (*seq.begin()).rest.return_type();
461 }
462    
463 unsigned add::return_type_tinfo(void) const
464 {
465     if (seq.size()==0) {
466         return tinfo_key;
467     }
468     return (*seq.begin()).rest.return_type_tinfo();
469 }
470
471 ex add::thisexpairseq(epvector const & v, ex const & oc) const
472 {
473     return (new add(v,oc))->setflag(status_flags::dynallocated);
474 }
475
476 ex add::thisexpairseq(epvector * vp, ex const & oc) const
477 {
478     return (new add(vp,oc))->setflag(status_flags::dynallocated);
479 }
480
481 /*
482 expair add::split_ex_to_pair(ex const & e) const
483 {
484     if (is_ex_exactly_of_type(e,mul)) {
485         mul const & mulref=ex_to_mul(e);
486         GINAC_ASSERT(mulref.seq.size()>1);
487         ex const & lastfactor_rest=(*(mulref.seq.end()-1)).rest;
488         ex const & lastfactor_coeff=(*(mulref.seq.end()-1)).coeff;
489         if (is_ex_exactly_of_type(lastfactor_rest,numeric) &&
490             ex_to_numeric(lastfactor_coeff).is_equal(numONE())) {
491             epvector s=mulref.seq;
492             //s.pop_back();
493             //return expair((new mul(s,1))->setflag(status_flags::dynallocated),
494             //              lastfactor);
495             mul * mulp=static_cast<mul *>(mulref.duplicate());
496 #ifdef EXPAIRSEQ_USE_HASHTAB
497             mulp->remove_hashtab_entry(mulp->seq.end()-1);
498 #endif // def EXPAIRSEQ_USE_HASHTAB
499             mulp->seq.pop_back();
500 #ifdef EXPAIRSEQ_USE_HASHTAB
501             mulp->shrink_hashtab();
502 #endif // def EXPAIRSEQ_USE_HASHTAB
503             mulp->clearflag(status_flags::evaluated);
504             mulp->clearflag(status_flags::hash_calculated);
505             return expair(mulp->setflag(status_flags::dynallocated),lastfactor_rest);
506         }
507     }
508     return expair(e,exONE());
509 }
510 */
511
512 expair add::split_ex_to_pair(ex const & e) const
513 {
514     if (is_ex_exactly_of_type(e,mul)) {
515         mul const & mulref=ex_to_mul(e);
516         ex numfactor=mulref.overall_coeff;
517         // mul * mulcopyp=static_cast<mul *>(mulref.duplicate());
518         mul * mulcopyp=new mul(mulref);
519         mulcopyp->overall_coeff=exONE();
520         mulcopyp->clearflag(status_flags::evaluated);
521         mulcopyp->clearflag(status_flags::hash_calculated);
522         return expair(mulcopyp->setflag(status_flags::dynallocated),numfactor);
523     }
524     return expair(e,exONE());
525 }
526
527 /*
528 expair add::combine_ex_with_coeff_to_pair(ex const & e,
529                                           ex const & c) const
530 {
531     GINAC_ASSERT(is_ex_exactly_of_type(c,numeric));
532     if (is_ex_exactly_of_type(e,mul)) {
533         mul const & mulref=ex_to_mul(e);
534         GINAC_ASSERT(mulref.seq.size()>1);
535         ex const & lastfactor_rest=(*(mulref.seq.end()-1)).rest;
536         ex const & lastfactor_coeff=(*(mulref.seq.end()-1)).coeff;
537         if (is_ex_exactly_of_type(lastfactor_rest,numeric) &&
538             ex_to_numeric(lastfactor_coeff).is_equal(numONE())) {
539             //epvector s=mulref.seq;
540             //s.pop_back();
541             //return expair((new mul(s,1))->setflag(status_flags::dynallocated),
542             //              ex_to_numeric(lastfactor).mul_dyn(ex_to_numeric(c)));
543             mul * mulp=static_cast<mul *>(mulref.duplicate());
544 #ifdef EXPAIRSEQ_USE_HASHTAB
545             mulp->remove_hashtab_entry(mulp->seq.end()-1);
546 #endif // def EXPAIRSEQ_USE_HASHTAB
547             mulp->seq.pop_back();
548 #ifdef EXPAIRSEQ_USE_HASHTAB
549             mulp->shrink_hashtab();
550 #endif // def EXPAIRSEQ_USE_HASHTAB
551             mulp->clearflag(status_flags::evaluated);
552             mulp->clearflag(status_flags::hash_calculated);
553             if (are_ex_trivially_equal(c,exONE())) {
554                 return expair(mulp->setflag(status_flags::dynallocated),lastfactor_rest);
555             } else if (are_ex_trivially_equal(lastfactor_rest,exONE())) {
556                 return expair(mulp->setflag(status_flags::dynallocated),c);
557             }                
558             return expair(mulp->setflag(status_flags::dynallocated),
559                           ex_to_numeric(lastfactor_rest).mul_dyn(ex_to_numeric(c)));
560         }
561     }
562     return expair(e,c);
563 }
564 */
565
566 expair add::combine_ex_with_coeff_to_pair(ex const & e,
567                                           ex const & c) const
568 {
569     GINAC_ASSERT(is_ex_exactly_of_type(c,numeric));
570     if (is_ex_exactly_of_type(e,mul)) {
571         mul const & mulref=ex_to_mul(e);
572         ex numfactor=mulref.overall_coeff;
573         //mul * mulcopyp=static_cast<mul *>(mulref.duplicate());
574         mul * mulcopyp=new mul(mulref);
575         mulcopyp->overall_coeff=exONE();
576         mulcopyp->clearflag(status_flags::evaluated);
577         mulcopyp->clearflag(status_flags::hash_calculated);
578         if (are_ex_trivially_equal(c,exONE())) {
579             return expair(mulcopyp->setflag(status_flags::dynallocated),numfactor);
580         } else if (are_ex_trivially_equal(numfactor,exONE())) {
581             return expair(mulcopyp->setflag(status_flags::dynallocated),c);
582         }
583         return expair(mulcopyp->setflag(status_flags::dynallocated),
584                           ex_to_numeric(numfactor).mul_dyn(ex_to_numeric(c)));
585     } else if (is_ex_exactly_of_type(e,numeric)) {
586         if (are_ex_trivially_equal(c,exONE())) {
587             return expair(e,exONE());
588         }
589         return expair(ex_to_numeric(e).mul_dyn(ex_to_numeric(c)),exONE());
590     }
591     return expair(e,c);
592 }
593     
594 expair add::combine_pair_with_coeff_to_pair(expair const & p,
595                                             ex const & c) const
596 {
597     GINAC_ASSERT(is_ex_exactly_of_type(p.coeff,numeric));
598     GINAC_ASSERT(is_ex_exactly_of_type(c,numeric));
599
600     if (is_ex_exactly_of_type(p.rest,numeric)) {
601         GINAC_ASSERT(ex_to_numeric(p.coeff).is_equal(numONE())); // should be normalized
602         return expair(ex_to_numeric(p.rest).mul_dyn(ex_to_numeric(c)),exONE());
603     }
604
605     return expair(p.rest,ex_to_numeric(p.coeff).mul_dyn(ex_to_numeric(c)));
606 }
607     
608 ex add::recombine_pair_to_ex(expair const & p) const
609 {
610     //if (p.coeff.compare(exONE())==0) {
611     //if (are_ex_trivially_equal(p.coeff,exONE())) {
612     if (ex_to_numeric(p.coeff).is_equal(numONE())) {
613         return p.rest;
614     } else {
615         return p.rest*p.coeff;
616     }
617 }
618
619 ex add::expand(unsigned options) const
620 {
621     epvector * vp=expandchildren(options);
622     if (vp==0) {
623         return *this;
624     }
625     return (new add(vp,overall_coeff))->setflag(status_flags::expanded    |
626                                                 status_flags::dynallocated );
627 }
628
629 //////////
630 // new virtual functions which can be overridden by derived classes
631 //////////
632
633 // none
634
635 //////////
636 // non-virtual functions in this class
637 //////////
638
639 // none
640
641 //////////
642 // static member variables
643 //////////
644
645 // protected
646
647 unsigned add::precedence=40;
648
649 //////////
650 // global constants
651 //////////
652
653 const add some_add;
654 type_info const & typeid_add=typeid(some_add);
655
656 } // namespace GiNaC