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