c899d19f1659c2641d3b4048108bf5b90019c5f7
[ginac.git] / ginac / basic.cpp
1 /** @file basic.cpp
2  *
3  *  Implementation of GiNaC's ABC. */
4
5 /*
6  *  GiNaC Copyright (C) 1999-2001 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 #ifdef DO_GINAC_ASSERT
26 #  include <typeinfo>
27 #endif
28
29 #include "basic.h"
30 #include "ex.h"
31 #include "numeric.h"
32 #include "power.h"
33 #include "symbol.h"
34 #include "lst.h"
35 #include "ncmul.h"
36 #include "relational.h"
37 #include "print.h"
38 #include "archive.h"
39 #include "utils.h"
40 #include "debugmsg.h"
41
42 namespace GiNaC {
43
44 GINAC_IMPLEMENT_REGISTERED_CLASS_NO_CTORS(basic, void)
45
46 //////////
47 // default ctor, dtor, copy ctor assignment operator and helpers
48 //////////
49
50 // public
51
52 basic::basic(const basic & other) : tinfo_key(TINFO_basic), flags(0), refcount(0)
53 {
54         debugmsg("basic copy ctor", LOGLEVEL_CONSTRUCT);
55         copy(other);
56 }
57
58 const basic & basic::operator=(const basic & other)
59 {
60         debugmsg("basic operator=", LOGLEVEL_ASSIGNMENT);
61         if (this != &other) {
62                 destroy(true);
63                 copy(other);
64         }
65         return *this;
66 }
67
68 // protected
69
70 // none (all conditionally inlined)
71
72 //////////
73 // other ctors
74 //////////
75
76 // none (all conditionally inlined)
77
78 //////////
79 // archiving
80 //////////
81
82 /** Construct object from archive_node. */
83 basic::basic(const archive_node &n, const lst &sym_lst) : flags(0), refcount(0)
84 {
85         debugmsg("basic ctor from archive_node", LOGLEVEL_CONSTRUCT);
86
87         // Reconstruct tinfo_key from class name
88         std::string class_name;
89         if (n.find_string("class", class_name))
90                 tinfo_key = find_tinfo_key(class_name);
91         else
92                 throw (std::runtime_error("archive node contains no class name"));
93 }
94
95 /** Unarchive the object. */
96 DEFAULT_UNARCHIVE(basic)
97
98 /** Archive the object. */
99 void basic::archive(archive_node &n) const
100 {
101         n.add_string("class", class_name());
102 }
103
104 //////////
105 // functions overriding virtual functions from bases classes
106 //////////
107
108 // none
109
110 //////////
111 // new virtual functions which can be overridden by derived classes
112 //////////
113
114 // public
115
116 /** Output to stream.
117  *  @param c print context object that describes the output formatting
118  *  @param level value that is used to identify the precedence or indentation
119  *               level for placing parentheses and formatting */
120 void basic::print(const print_context & c, unsigned level) const
121 {
122         debugmsg("basic print", LOGLEVEL_PRINT);
123
124         if (is_of_type(c, print_tree)) {
125
126                 c.s << std::string(level, ' ') << class_name()
127                     << std::hex << ", hash=0x" << hashvalue << ", flags=0x" << flags << std::dec
128                     << ", nops=" << nops()
129                     << std::endl;
130                 for (unsigned i=0; i<nops(); ++i)
131                         op(i).print(c, level + static_cast<const print_tree &>(c).delta_indent);
132
133         } else
134                 c.s << "[" << class_name() << " object]";
135 }
136
137 /** Little wrapper arount print to be called within a debugger.
138  *  This is needed because you cannot call foo.print(cout) from within the
139  *  debugger because it might not know what cout is.  This method can be
140  *  invoked with no argument and it will simply print to stdout.
141  *
142  *  @see basic::print */
143 void basic::dbgprint(void) const
144 {
145         this->print(std::cerr);
146         std::cerr << std::endl;
147 }
148
149 /** Little wrapper arount printtree to be called within a debugger.
150  *
151  *  @see basic::dbgprint
152  *  @see basic::printtree */
153 void basic::dbgprinttree(void) const
154 {
155         this->print(print_tree(std::cerr));
156 }
157
158 /** Return relative operator precedence (for parenthizing output). */
159 unsigned basic::precedence(void) const
160 {
161         return 70;
162 }
163
164 /** Create a new copy of this on the heap.  One can think of this as simulating
165  *  a virtual copy constructor which is needed for instance by the refcounted
166  *  construction of an ex from a basic. */
167 basic * basic::duplicate() const
168 {
169         debugmsg("basic duplicate",LOGLEVEL_DUPLICATE);
170         return new basic(*this);
171 }
172
173 /** Information about the object.
174  *
175  *  @see class info_flags */
176 bool basic::info(unsigned inf) const
177 {
178         // all possible properties are false for basic objects
179         return false;
180 }
181
182 /** Number of operands/members. */
183 unsigned basic::nops() const
184 {
185         // iterating from 0 to nops() on atomic objects should be an empty loop,
186         // and accessing their elements is a range error.  Container objects should
187         // override this.
188         return 0;
189 }
190
191 /** Return operand/member at position i. */
192 ex basic::op(int i) const
193 {
194         return (const_cast<basic *>(this))->let_op(i);
195 }
196
197 /** Return modifyable operand/member at position i. */
198 ex & basic::let_op(int i)
199 {
200         throw(std::out_of_range("op() out of range"));
201 }
202
203 ex basic::operator[](const ex & index) const
204 {
205         if (is_exactly_of_type(*index.bp,numeric))
206                 return op(static_cast<const numeric &>(*index.bp).to_int());
207         
208         throw(std::invalid_argument("non-numeric indices not supported by this type"));
209 }
210
211 ex basic::operator[](int i) const
212 {
213         return op(i);
214 }
215
216 /** Search ocurrences.  An object  'has' an expression if it is the expression
217  *  itself or one of the children 'has' it.  As a consequence (according to
218  *  the definition of children) given e=x+y+z, e.has(x) is true but e.has(x+y)
219  *  is false. */
220 bool basic::has(const ex & other) const
221 {
222         GINAC_ASSERT(other.bp!=0);
223         lst repl_lst;
224         if (match(*other.bp, repl_lst)) return true;
225         if (nops()>0) {
226                 for (unsigned i=0; i<nops(); i++)
227                         if (op(i).has(other))
228                                 return true;
229         }
230         
231         return false;
232 }
233
234 /** Return degree of highest power in object s. */
235 int basic::degree(const ex & s) const
236 {
237         return 0;
238 }
239
240 /** Return degree of lowest power in object s. */
241 int basic::ldegree(const ex & s) const
242 {
243         return 0;
244 }
245
246 /** Return coefficient of degree n in object s. */
247 ex basic::coeff(const ex & s, int n) const
248 {
249         return n==0 ? *this : _ex0();
250 }
251
252 /** Sort expression in terms of powers of some object(s).
253  *  @param s object(s) to sort in
254  *  @param distributed recursive or distributed form (only used when s is a list) */
255 ex basic::collect(const ex & s, bool distributed) const
256 {
257         ex x;
258         if (is_ex_of_type(s, lst)) {
259
260                 // List of objects specified
261                 if (s.nops() == 1)
262                         return collect(s.op(0));
263
264                 else if (distributed) {
265
266                         // Get lower/upper degree of all symbols in list
267                         int num = s.nops();
268                         struct sym_info {
269                                 ex sym;
270                                 int ldeg, deg;
271                                 int cnt;  // current degree, 'counter'
272                                 ex coeff; // coefficient for degree 'cnt'
273                         };
274                         sym_info *si = new sym_info[num];
275                         ex c = *this;
276                         for (int i=0; i<num; i++) {
277                                 si[i].sym = s.op(i);
278                                 si[i].ldeg = si[i].cnt = this->ldegree(si[i].sym);
279                                 si[i].deg = this->degree(si[i].sym);
280                                 c = si[i].coeff = c.coeff(si[i].sym, si[i].cnt);
281                         }
282
283                         while (true) {
284
285                                 // Calculate coeff*x1^c1*...*xn^cn
286                                 ex y = _ex1();
287                                 for (int i=0; i<num; i++) {
288                                         int cnt = si[i].cnt;
289                                         y *= power(si[i].sym, cnt);
290                                 }
291                                 x += y * si[num - 1].coeff;
292
293                                 // Increment counters
294                                 int n = num - 1;
295                                 while (true) {
296                                         si[n].cnt++;
297                                         if (si[n].cnt <= si[n].deg) {
298                                                 // Update coefficients
299                                                 ex c;
300                                                 if (n == 0)
301                                                         c = *this;
302                                                 else
303                                                         c = si[n - 1].coeff;
304                                                 for (int i=n; i<num; i++)
305                                                         c = si[i].coeff = c.coeff(si[i].sym, si[i].cnt);
306                                                 break;
307                                         }
308                                         if (n == 0)
309                                                 goto done;
310                                         si[n].cnt = si[n].ldeg;
311                                         n--;
312                                 }
313                         }
314
315 done:           delete[] si;
316
317                 } else {
318
319                         // Recursive form
320                         x = *this;
321                         for (int n=s.nops()-1; n>=0; n--)
322                                 x = x.collect(s[n]);
323                 }
324
325         } else {
326
327                 // Only one object specified
328                 for (int n=this->ldegree(s); n<=this->degree(s); ++n)
329                         x += this->coeff(s,n)*power(s,n);
330         }
331         
332         // correct for lost fractional arguments and return
333         return x + (*this - x).expand();
334 }
335
336 /** Perform automatic non-interruptive symbolic evaluation on expression. */
337 ex basic::eval(int level) const
338 {
339         // There is nothing to do for basic objects:
340         return this->hold();
341 }
342
343 /** Evaluate object numerically. */
344 ex basic::evalf(int level) const
345 {
346         // There is nothing to do for basic objects:
347         return *this;
348 }
349
350 /** Perform automatic symbolic evaluations on indexed expression that
351  *  contains this object as the base expression. */
352 ex basic::eval_indexed(const basic & i) const
353  // this function can't take a "const ex & i" because that would result
354  // in an infinite eval() loop
355 {
356         // There is nothing to do for basic objects
357         return i.hold();
358 }
359
360 /** Add two indexed expressions. They are guaranteed to be of class indexed
361  *  (or a subclass) and their indices are compatible. This function is used
362  *  internally by simplify_indexed().
363  *
364  *  @param self First indexed expression; it's base object is *this
365  *  @param other Second indexed expression
366  *  @return sum of self and other 
367  *  @see ex::simplify_indexed() */
368 ex basic::add_indexed(const ex & self, const ex & other) const
369 {
370         return self + other;
371 }
372
373 /** Multiply an indexed expression with a scalar. This function is used
374  *  internally by simplify_indexed().
375  *
376  *  @param self Indexed expression; it's base object is *this
377  *  @param other Numeric value
378  *  @return product of self and other
379  *  @see ex::simplify_indexed() */
380 ex basic::scalar_mul_indexed(const ex & self, const numeric & other) const
381 {
382         return self * other;
383 }
384
385 /** Try to contract two indexed expressions that appear in the same product. 
386  *  If a contraction exists, the function overwrites one or both of the
387  *  expressions and returns true. Otherwise it returns false. It is
388  *  guaranteed that both expressions are of class indexed (or a subclass)
389  *  and that at least one dummy index has been found. This functions is
390  *  used internally by simplify_indexed().
391  *
392  *  @param self Pointer to first indexed expression; it's base object is *this
393  *  @param other Pointer to second indexed expression
394  *  @param v The complete vector of factors
395  *  @return true if the contraction was successful, false otherwise
396  *  @see ex::simplify_indexed() */
397 bool basic::contract_with(exvector::iterator self, exvector::iterator other, exvector & v) const
398 {
399         // Do nothing
400         return false;
401 }
402
403 /** Check whether the expression matches a given pattern. For every wildcard
404  *  object in the pattern, an expression of the form "wildcard == matching_expression"
405  *  is added to repl_lst. */
406 bool basic::match(const ex & pattern, lst & repl_lst) const
407 {
408 //clog << "match " << *this << " with " << pattern << ", repl_lst = " << repl_lst << endl;
409         if (is_ex_exactly_of_type(pattern, wildcard)) {
410
411                 // Wildcard matches anything, but check whether we already have found
412                 // a match for that wildcard first (if so, it the earlier match must
413                 // be the same expression)
414                 for (unsigned i=0; i<repl_lst.nops(); i++) {
415                         if (repl_lst.op(i).op(0).is_equal(pattern))
416                                 return is_equal(*repl_lst.op(i).op(1).bp);
417                 }
418                 repl_lst.append(pattern == *this);
419                 return true;
420
421         } else {
422
423                 // Expression must be of the same type as the pattern
424                 if (tinfo() != pattern.bp->tinfo())
425                         return false;
426
427                 // Number of subexpressions must match
428                 if (nops() != pattern.nops())
429                         return false;
430
431                 // No subexpressions? Then just compare the objects (there can't be
432                 // wildcards in the pattern)
433                 if (nops() == 0)
434                         return is_equal(*pattern.bp);
435
436                 // Otherwise the subexpressions must match one-to-one
437                 for (unsigned i=0; i<nops(); i++)
438                         if (!op(i).match(pattern.op(i), repl_lst))
439                                 return false;
440
441                 // Looks similar enough, match found
442                 return true;
443         }
444 }
445
446 /** Substitute a set of objects by arbitrary expressions. The ex returned
447  *  will already be evaluated. */
448 ex basic::subs(const lst & ls, const lst & lr, bool no_pattern) const
449 {
450         GINAC_ASSERT(ls.nops() == lr.nops());
451
452         if (no_pattern) {
453                 for (unsigned i=0; i<ls.nops(); i++) {
454                         if (is_equal(*ls.op(i).bp))
455                                 return lr.op(i);
456                 }
457         } else {
458                 for (unsigned i=0; i<ls.nops(); i++) {
459                         lst repl_lst;
460                         if (match(*ls.op(i).bp, repl_lst))
461                                 return lr.op(i).bp->subs(repl_lst, true); // avoid recursion when re-substituting the wildcards
462                 }
463         }
464
465         return *this;
466 }
467
468 /** Default interface of nth derivative ex::diff(s, n).  It should be called
469  *  instead of ::derivative(s) for first derivatives and for nth derivatives it
470  *  just recurses down.
471  *
472  *  @param s symbol to differentiate in
473  *  @param nth order of differentiation
474  *  @see ex::diff */
475 ex basic::diff(const symbol & s, unsigned nth) const
476 {
477         // trivial: zeroth derivative
478         if (nth==0)
479                 return ex(*this);
480         
481         // evaluate unevaluated *this before differentiating
482         if (!(flags & status_flags::evaluated))
483                 return ex(*this).diff(s, nth);
484         
485         ex ndiff = this->derivative(s);
486         while (!ndiff.is_zero() &&    // stop differentiating zeros
487                nth>1) {
488                 ndiff = ndiff.diff(s);
489                 --nth;
490         }
491         return ndiff;
492 }
493
494 /** Return a vector containing the free indices of an expression. */
495 exvector basic::get_free_indices(void) const
496 {
497         return exvector(); // return an empty exvector
498 }
499
500 ex basic::simplify_ncmul(const exvector & v) const
501 {
502         return simplified_ncmul(v);
503 }
504
505 // protected
506
507 /** Default implementation of ex::diff(). It simply throws an error message.
508  *
509  *  @exception logic_error (differentiation not supported by this type)
510  *  @see ex::diff */
511 ex basic::derivative(const symbol & s) const
512 {
513         throw(std::logic_error("differentiation not supported by this type"));
514 }
515
516 /** Returns order relation between two objects of same type.  This needs to be
517  *  implemented by each class. It may never return anything else than 0,
518  *  signalling equality, or +1 and -1 signalling inequality and determining
519  *  the canonical ordering.  (Perl hackers will wonder why C++ doesn't feature
520  *  the spaceship operator <=> for denoting just this.) */
521 int basic::compare_same_type(const basic & other) const
522 {
523         return compare_pointers(this, &other);
524 }
525
526 /** Returns true if two objects of same type are equal.  Normally needs
527  *  not be reimplemented as long as it wasn't overwritten by some parent
528  *  class, since it just calls compare_same_type().  The reason why this
529  *  function exists is that sometimes it is easier to determine equality
530  *  than an order relation and then it can be overridden. */
531 bool basic::is_equal_same_type(const basic & other) const
532 {
533         return this->compare_same_type(other)==0;
534 }
535
536 unsigned basic::return_type(void) const
537 {
538         return return_types::commutative;
539 }
540
541 unsigned basic::return_type_tinfo(void) const
542 {
543         return tinfo();
544 }
545
546 /** Compute the hash value of an object and if it makes sense to store it in
547  *  the objects status_flags, do so.  The method inherited from class basic
548  *  computes a hash value based on the type and hash values of possible
549  *  members.  For this reason it is well suited for container classes but
550  *  atomic classes should override this implementation because otherwise they
551  *  would all end up with the same hashvalue. */
552 unsigned basic::calchash(void) const
553 {
554         unsigned v = golden_ratio_hash(tinfo());
555         for (unsigned i=0; i<nops(); i++) {
556                 v = rotate_left_31(v);
557                 v ^= (const_cast<basic *>(this))->op(i).gethash();
558         }
559         
560         // mask out numeric hashes:
561         v &= 0x7FFFFFFFU;
562         
563         // store calculated hash value only if object is already evaluated
564         if (flags & status_flags::evaluated) {
565                 setflag(status_flags::hash_calculated);
566                 hashvalue = v;
567         }
568
569         return v;
570 }
571
572 /** Expand expression, i.e. multiply it out and return the result as a new
573  *  expression. */
574 ex basic::expand(unsigned options) const
575 {
576         return this->setflag(status_flags::expanded);
577 }
578
579
580 //////////
581 // non-virtual functions in this class
582 //////////
583
584 // public
585
586 /** Substitute objects in an expression (syntactic substitution) and return
587  *  the result as a new expression.  There are two valid types of
588  *  replacement arguments: 1) a relational like object==ex and 2) a list of
589  *  relationals lst(object1==ex1,object2==ex2,...), which is converted to
590  *  subs(lst(object1,object2,...),lst(ex1,ex2,...)). */
591 ex basic::subs(const ex & e, bool no_pattern) const
592 {
593         if (e.info(info_flags::relation_equal)) {
594                 return subs(lst(e), no_pattern);
595         }
596         if (!e.info(info_flags::list)) {
597                 throw(std::invalid_argument("basic::subs(ex): argument must be a list"));
598         }
599         lst ls;
600         lst lr;
601         for (unsigned i=0; i<e.nops(); i++) {
602                 ex r = e.op(i);
603                 if (!r.info(info_flags::relation_equal)) {
604                         throw(std::invalid_argument("basic::subs(ex): argument must be a list of equations"));
605                 }
606                 ls.append(r.op(0));
607                 lr.append(r.op(1));
608         }
609         return subs(ls, lr, no_pattern);
610 }
611
612 /** Compare objects to establish canonical ordering.
613  *  All compare functions return: -1 for *this less than other, 0 equal,
614  *  1 greater. */
615 int basic::compare(const basic & other) const
616 {
617         unsigned hash_this = gethash();
618         unsigned hash_other = other.gethash();
619         
620         if (hash_this<hash_other) return -1;
621         if (hash_this>hash_other) return 1;
622         
623         unsigned typeid_this = tinfo();
624         unsigned typeid_other = other.tinfo();
625         
626         if (typeid_this<typeid_other) {
627 //              std::cout << "hash collision, different types: " 
628 //                        << *this << " and " << other << std::endl;
629 //              this->print(print_tree(std::cout));
630 //              std::cout << " and ";
631 //              other.print(print_tree(std::cout));
632 //              std::cout << std::endl;
633                 return -1;
634         }
635         if (typeid_this>typeid_other) {
636 //              std::cout << "hash collision, different types: " 
637 //                        << *this << " and " << other << std::endl;
638 //              this->print(print_tree(std::cout));
639 //              std::cout << " and ";
640 //              other.print(print_tree(std::cout));
641 //              std::cout << std::endl;
642                 return 1;
643         }
644         
645         GINAC_ASSERT(typeid(*this)==typeid(other));
646         
647 //      int cmpval = compare_same_type(other);
648 //      if ((cmpval!=0) && (hash_this<0x80000000U)) {
649 //              std::cout << "hash collision, same type: " 
650 //                        << *this << " and " << other << std::endl;
651 //              this->print(print_tree(std::cout));
652 //              std::cout << " and ";
653 //              other.print(print_tree(std::cout));
654 //              std::cout << std::endl;
655 //      }
656 //      return cmpval;
657         
658         return compare_same_type(other);
659 }
660
661 /** Test for equality.
662  *  This is only a quick test, meaning objects should be in the same domain.
663  *  You might have to .expand(), .normal() objects first, depending on the
664  *  domain of your computation, to get a more reliable answer.
665  *
666  *  @see is_equal_same_type */
667 bool basic::is_equal(const basic & other) const
668 {
669         if (this->gethash()!=other.gethash())
670                 return false;
671         if (this->tinfo()!=other.tinfo())
672                 return false;
673         
674         GINAC_ASSERT(typeid(*this)==typeid(other));
675         
676         return this->is_equal_same_type(other);
677 }
678
679 // protected
680
681 /** Stop further evaluation.
682  *
683  *  @see basic::eval */
684 const basic & basic::hold(void) const
685 {
686         return this->setflag(status_flags::evaluated);
687 }
688
689 /** Ensure the object may be modified without hurting others, throws if this
690  *  is not the case. */
691 void basic::ensure_if_modifiable(void) const
692 {
693         if (this->refcount>1)
694                 throw(std::runtime_error("cannot modify multiply referenced object"));
695 }
696
697 //////////
698 // global variables
699 //////////
700
701 int max_recursion_level = 1024;
702
703 } // namespace GiNaC