]> www.ginac.de Git - ginac.git/blob - ginac/basic.cpp
725f28a2e6b0ea40965ebbc45d50e6ca9d5e7a98
[ginac.git] / ginac / basic.cpp
1 /** @file basic.cpp
2  *
3  *  Implementation of GiNaC's ABC. */
4
5 /*
6  *  GiNaC Copyright (C) 1999-2003 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 "operators.h"
38 #include "wildcard.h"
39 #include "archive.h"
40 #include "utils.h"
41
42 namespace GiNaC {
43
44 GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(basic, void,
45   print_func<print_context>(&basic::do_print).
46   print_func<print_tree>(&basic::do_print_tree).
47   print_func<print_python_repr>(&basic::do_print_python_repr))
48
49 //////////
50 // default constructor, destructor, copy constructor and assignment operator
51 //////////
52
53 // public
54
55 /** basic copy constructor: implicitly assumes that the other class is of
56  *  the exact same type (as it's used by duplicate()), so it can copy the
57  *  tinfo_key and the hash value. */
58 basic::basic(const basic & other) : tinfo_key(other.tinfo_key), flags(other.flags & ~status_flags::dynallocated), hashvalue(other.hashvalue), refcount(0)
59 {
60         GINAC_ASSERT(typeid(*this) == typeid(other));
61 }
62
63 /** basic assignment operator: the other object might be of a derived class. */
64 const basic & basic::operator=(const basic & other)
65 {
66         unsigned fl = other.flags & ~status_flags::dynallocated;
67         if (tinfo_key != other.tinfo_key) {
68                 // The other object is of a derived class, so clear the flags as they
69                 // might no longer apply (especially hash_calculated). Oh, and don't
70                 // copy the tinfo_key: it is already set correctly for this object.
71                 flags = 0;
72         } else {
73                 // The objects are of the exact same class, so copy the hash value.
74                 hashvalue = other.hashvalue;
75         }
76         flags = fl;
77         refcount = 0;
78         return *this;
79 }
80
81 // protected
82
83 // none (all inlined)
84
85 //////////
86 // other constructors
87 //////////
88
89 // none (all inlined)
90
91 //////////
92 // archiving
93 //////////
94
95 /** Construct object from archive_node. */
96 basic::basic(const archive_node &n, lst &sym_lst) : flags(0), refcount(0)
97 {
98         // Reconstruct tinfo_key from class name
99         std::string class_name;
100         if (n.find_string("class", class_name))
101                 tinfo_key = find_tinfo_key(class_name);
102         else
103                 throw (std::runtime_error("archive node contains no class name"));
104 }
105
106 /** Unarchive the object. */
107 DEFAULT_UNARCHIVE(basic)
108
109 /** Archive the object. */
110 void basic::archive(archive_node &n) const
111 {
112         n.add_string("class", class_name());
113 }
114
115 //////////
116 // new virtual functions which can be overridden by derived classes
117 //////////
118
119 // public
120
121 /** Output to stream. This performs double dispatch on the dynamic type of
122  *  *this and the dynamic type of the supplied print context.
123  *  @param c print context object that describes the output formatting
124  *  @param level value that is used to identify the precedence or indentation
125  *               level for placing parentheses and formatting */
126 void basic::print(const print_context & c, unsigned level) const
127 {
128         print_dispatch(get_class_info(), c, level);
129 }
130
131 /** Like print(), but dispatch to the specified class. Can be used by
132  *  implementations of print methods to dispatch to the method of the
133  *  superclass.
134  *
135  *  @see basic::print */
136 void basic::print_dispatch(const registered_class_info & ri, const print_context & c, unsigned level) const
137 {
138         // Double dispatch on object type and print_context type
139         const registered_class_info * reg_info = &ri;
140         const print_context_class_info * pc_info = &c.get_class_info();
141
142 next_class:
143 std::clog << "searching class " << reg_info->options.get_name() << std::endl;
144         const std::vector<print_functor> & pdt = reg_info->options.get_print_dispatch_table();
145
146 next_context:
147 std::clog << "searching context " << pc_info->options.get_name() << ", ID " << pc_info->options.get_id() << std::endl;
148         unsigned id = pc_info->options.get_id();
149         if (id >= pdt.size() || !(pdt[id].is_valid())) {
150
151                 // Method not found, try parent print_context class
152                 const print_context_class_info * parent_pc_info = pc_info->get_parent();
153                 if (parent_pc_info) {
154                         pc_info = parent_pc_info;
155                         goto next_context;
156                 }
157
158                 // Method still not found, try parent class
159                 const registered_class_info * parent_reg_info = reg_info->get_parent();
160                 if (parent_reg_info) {
161                         reg_info = parent_reg_info;
162                         pc_info = &c.get_class_info();
163                         goto next_class;
164                 }
165
166                 // Method still not found. This shouldn't happen because basic (the
167                 // base class of the algebraic hierarchy) registers a method for
168                 // print_context (the base class of the print context hierarchy),
169                 // so if we end up here, there's something wrong with the class
170                 // registry.
171                 throw (std::runtime_error(std::string("basic::print(): method for ") + class_name() + "/" + c.class_name() + " not found"));
172
173         } else {
174
175                 // Call method
176 std::clog << "method found, calling" << std::endl;
177 std::clog << " this = " << class_name() << ", context = " << c.class_name() << std::endl;
178                 pdt[id](*this, c, level);
179         }
180 }
181
182 /** Default output to stream. */
183 void basic::do_print(const print_context & c, unsigned level) const
184 {
185         c.s << "[" << class_name() << " object]";
186 }
187
188 /** Tree output to stream. */
189 void basic::do_print_tree(const print_tree & c, unsigned level) const
190 {
191         c.s << std::string(level, ' ') << class_name()
192             << std::hex << ", hash=0x" << hashvalue << ", flags=0x" << flags << std::dec
193             << ", nops=" << nops()
194             << std::endl;
195         for (size_t i=0; i<nops(); ++i)
196                 op(i).print(c, level + c.delta_indent);
197 }
198
199 /** Python parsable output to stream. */
200 void basic::do_print_python_repr(const print_python_repr & c, unsigned level) const
201 {
202         c.s << class_name() << "()";
203 }
204
205 /** Little wrapper around print to be called within a debugger.
206  *  This is needed because you cannot call foo.print(cout) from within the
207  *  debugger because it might not know what cout is.  This method can be
208  *  invoked with no argument and it will simply print to stdout.
209  *
210  *  @see basic::print
211  *  @see basic::dbgprinttree */
212 void basic::dbgprint() const
213 {
214         this->print(std::cerr);
215         std::cerr << std::endl;
216 }
217
218 /** Little wrapper around printtree to be called within a debugger.
219  *
220  *  @see basic::dbgprint */
221 void basic::dbgprinttree() const
222 {
223         this->print(print_tree(std::cerr));
224 }
225
226 /** Return relative operator precedence (for parenthezing output). */
227 unsigned basic::precedence() const
228 {
229         return 70;
230 }
231
232 /** Information about the object.
233  *
234  *  @see class info_flags */
235 bool basic::info(unsigned inf) const
236 {
237         // all possible properties are false for basic objects
238         return false;
239 }
240
241 /** Number of operands/members. */
242 size_t basic::nops() const
243 {
244         // iterating from 0 to nops() on atomic objects should be an empty loop,
245         // and accessing their elements is a range error.  Container objects should
246         // override this.
247         return 0;
248 }
249
250 /** Return operand/member at position i. */
251 ex basic::op(size_t i) const
252 {
253         throw(std::range_error(std::string("basic::op(): ") + class_name() + std::string(" has no operands")));
254 }
255
256 /** Return modifyable operand/member at position i. */
257 ex & basic::let_op(size_t i)
258 {
259         ensure_if_modifiable();
260         throw(std::range_error(std::string("basic::let_op(): ") + class_name() + std::string(" has no operands")));
261 }
262
263 ex basic::operator[](const ex & index) const
264 {
265         if (is_exactly_a<numeric>(index))
266                 return op(static_cast<size_t>(ex_to<numeric>(index).to_int()));
267
268         throw(std::invalid_argument(std::string("non-numeric indices not supported by ") + class_name()));
269 }
270
271 ex basic::operator[](size_t i) const
272 {
273         return op(i);
274 }
275
276 ex & basic::operator[](const ex & index)
277 {
278         if (is_exactly_a<numeric>(index))
279                 return let_op(ex_to<numeric>(index).to_int());
280
281         throw(std::invalid_argument(std::string("non-numeric indices not supported by ") + class_name()));
282 }
283
284 ex & basic::operator[](size_t i)
285 {
286         return let_op(i);
287 }
288
289 /** Test for occurrence of a pattern.  An object 'has' a pattern if it matches
290  *  the pattern itself or one of the children 'has' it.  As a consequence
291  *  (according to the definition of children) given e=x+y+z, e.has(x) is true
292  *  but e.has(x+y) is false. */
293 bool basic::has(const ex & pattern) const
294 {
295         lst repl_lst;
296         if (match(pattern, repl_lst))
297                 return true;
298         for (size_t i=0; i<nops(); i++)
299                 if (op(i).has(pattern))
300                         return true;
301         
302         return false;
303 }
304
305 /** Construct new expression by applying the specified function to all
306  *  sub-expressions (one level only, not recursively). */
307 ex basic::map(map_function & f) const
308 {
309         size_t num = nops();
310         if (num == 0)
311                 return *this;
312
313         basic *copy = duplicate();
314         copy->setflag(status_flags::dynallocated);
315         copy->clearflag(status_flags::hash_calculated | status_flags::expanded);
316         for (size_t i=0; i<num; i++)
317                 copy->let_op(i) = f(copy->op(i));
318         return *copy;
319 }
320
321 /** Return degree of highest power in object s. */
322 int basic::degree(const ex & s) const
323 {
324         return is_equal(ex_to<basic>(s)) ? 1 : 0;
325 }
326
327 /** Return degree of lowest power in object s. */
328 int basic::ldegree(const ex & s) const
329 {
330         return is_equal(ex_to<basic>(s)) ? 1 : 0;
331 }
332
333 /** Return coefficient of degree n in object s. */
334 ex basic::coeff(const ex & s, int n) const
335 {
336         if (is_equal(ex_to<basic>(s)))
337                 return n==1 ? _ex1 : _ex0;
338         else
339                 return n==0 ? *this : _ex0;
340 }
341
342 /** Sort expanded expression in terms of powers of some object(s).
343  *  @param s object(s) to sort in
344  *  @param distributed recursive or distributed form (only used when s is a list) */
345 ex basic::collect(const ex & s, bool distributed) const
346 {
347         ex x;
348         if (is_a<lst>(s)) {
349
350                 // List of objects specified
351                 if (s.nops() == 0)
352                         return *this;
353                 if (s.nops() == 1)
354                         return collect(s.op(0));
355
356                 else if (distributed) {
357
358                         // Get lower/upper degree of all symbols in list
359                         size_t num = s.nops();
360                         struct sym_info {
361                                 ex sym;
362                                 int ldeg, deg;
363                                 int cnt;  // current degree, 'counter'
364                                 ex coeff; // coefficient for degree 'cnt'
365                         };
366                         sym_info *si = new sym_info[num];
367                         ex c = *this;
368                         for (size_t i=0; i<num; i++) {
369                                 si[i].sym = s.op(i);
370                                 si[i].ldeg = si[i].cnt = this->ldegree(si[i].sym);
371                                 si[i].deg = this->degree(si[i].sym);
372                                 c = si[i].coeff = c.coeff(si[i].sym, si[i].cnt);
373                         }
374
375                         while (true) {
376
377                                 // Calculate coeff*x1^c1*...*xn^cn
378                                 ex y = _ex1;
379                                 for (size_t i=0; i<num; i++) {
380                                         int cnt = si[i].cnt;
381                                         y *= power(si[i].sym, cnt);
382                                 }
383                                 x += y * si[num - 1].coeff;
384
385                                 // Increment counters
386                                 size_t n = num - 1;
387                                 while (true) {
388                                         ++si[n].cnt;
389                                         if (si[n].cnt <= si[n].deg) {
390                                                 // Update coefficients
391                                                 ex c;
392                                                 if (n == 0)
393                                                         c = *this;
394                                                 else
395                                                         c = si[n - 1].coeff;
396                                                 for (size_t i=n; i<num; i++)
397                                                         c = si[i].coeff = c.coeff(si[i].sym, si[i].cnt);
398                                                 break;
399                                         }
400                                         if (n == 0)
401                                                 goto done;
402                                         si[n].cnt = si[n].ldeg;
403                                         n--;
404                                 }
405                         }
406
407 done:           delete[] si;
408
409                 } else {
410
411                         // Recursive form
412                         x = *this;
413                         size_t n = s.nops() - 1;
414                         while (true) {
415                                 x = x.collect(s[n]);
416                                 if (n == 0)
417                                         break;
418                                 n--;
419                         }
420                 }
421
422         } else {
423
424                 // Only one object specified
425                 for (int n=this->ldegree(s); n<=this->degree(s); ++n)
426                         x += this->coeff(s,n)*power(s,n);
427         }
428         
429         // correct for lost fractional arguments and return
430         return x + (*this - x).expand();
431 }
432
433 /** Perform automatic non-interruptive term rewriting rules. */
434 ex basic::eval(int level) const
435 {
436         // There is nothing to do for basic objects:
437         return hold();
438 }
439
440 /** Function object to be applied by basic::evalf(). */
441 struct evalf_map_function : public map_function {
442         int level;
443         evalf_map_function(int l) : level(l) {}
444         ex operator()(const ex & e) { return evalf(e, level); }
445 };
446
447 /** Evaluate object numerically. */
448 ex basic::evalf(int level) const
449 {
450         if (nops() == 0)
451                 return *this;
452         else {
453                 if (level == 1)
454                         return *this;
455                 else if (level == -max_recursion_level)
456                         throw(std::runtime_error("max recursion level reached"));
457                 else {
458                         evalf_map_function map_evalf(level - 1);
459                         return map(map_evalf);
460                 }
461         }
462 }
463
464 /** Function object to be applied by basic::evalm(). */
465 struct evalm_map_function : public map_function {
466         ex operator()(const ex & e) { return evalm(e); }
467 } map_evalm;
468
469 /** Evaluate sums, products and integer powers of matrices. */
470 ex basic::evalm() const
471 {
472         if (nops() == 0)
473                 return *this;
474         else
475                 return map(map_evalm);
476 }
477
478 /** Perform automatic symbolic evaluations on indexed expression that
479  *  contains this object as the base expression. */
480 ex basic::eval_indexed(const basic & i) const
481  // this function can't take a "const ex & i" because that would result
482  // in an infinite eval() loop
483 {
484         // There is nothing to do for basic objects
485         return i.hold();
486 }
487
488 /** Add two indexed expressions. They are guaranteed to be of class indexed
489  *  (or a subclass) and their indices are compatible. This function is used
490  *  internally by simplify_indexed().
491  *
492  *  @param self First indexed expression; it's base object is *this
493  *  @param other Second indexed expression
494  *  @return sum of self and other 
495  *  @see ex::simplify_indexed() */
496 ex basic::add_indexed(const ex & self, const ex & other) const
497 {
498         return self + other;
499 }
500
501 /** Multiply an indexed expression with a scalar. This function is used
502  *  internally by simplify_indexed().
503  *
504  *  @param self Indexed expression; it's base object is *this
505  *  @param other Numeric value
506  *  @return product of self and other
507  *  @see ex::simplify_indexed() */
508 ex basic::scalar_mul_indexed(const ex & self, const numeric & other) const
509 {
510         return self * other;
511 }
512
513 /** Try to contract two indexed expressions that appear in the same product. 
514  *  If a contraction exists, the function overwrites one or both of the
515  *  expressions and returns true. Otherwise it returns false. It is
516  *  guaranteed that both expressions are of class indexed (or a subclass)
517  *  and that at least one dummy index has been found. This functions is
518  *  used internally by simplify_indexed().
519  *
520  *  @param self Pointer to first indexed expression; it's base object is *this
521  *  @param other Pointer to second indexed expression
522  *  @param v The complete vector of factors
523  *  @return true if the contraction was successful, false otherwise
524  *  @see ex::simplify_indexed() */
525 bool basic::contract_with(exvector::iterator self, exvector::iterator other, exvector & v) const
526 {
527         // Do nothing
528         return false;
529 }
530
531 /** Check whether the expression matches a given pattern. For every wildcard
532  *  object in the pattern, an expression of the form "wildcard == matching_expression"
533  *  is added to repl_lst. */
534 bool basic::match(const ex & pattern, lst & repl_lst) const
535 {
536 /*
537         Sweet sweet shapes, sweet sweet shapes,
538         That's the key thing, right right.
539         Feed feed face, feed feed shapes,
540         But who is the king tonight?
541         Who is the king tonight?
542         Pattern is the thing, the key thing-a-ling,
543         But who is the king of Pattern?
544         But who is the king, the king thing-a-ling,
545         Who is the king of Pattern?
546         Bog is the king, the king thing-a-ling,
547         Bog is the king of Pattern.
548         Ba bu-bu-bu-bu bu-bu-bu-bu-bu-bu bu-bu
549         Bog is the king of Pattern.
550 */
551
552         if (is_exactly_a<wildcard>(pattern)) {
553
554                 // Wildcard matches anything, but check whether we already have found
555                 // a match for that wildcard first (if so, it the earlier match must
556                 // be the same expression)
557                 for (lst::const_iterator it = repl_lst.begin(); it != repl_lst.end(); ++it) {
558                         if (it->op(0).is_equal(pattern))
559                                 return is_equal(ex_to<basic>(it->op(1)));
560                 }
561                 repl_lst.append(pattern == *this);
562                 return true;
563
564         } else {
565
566                 // Expression must be of the same type as the pattern
567                 if (tinfo() != ex_to<basic>(pattern).tinfo())
568                         return false;
569
570                 // Number of subexpressions must match
571                 if (nops() != pattern.nops())
572                         return false;
573
574                 // No subexpressions? Then just compare the objects (there can't be
575                 // wildcards in the pattern)
576                 if (nops() == 0)
577                         return is_equal_same_type(ex_to<basic>(pattern));
578
579                 // Check whether attributes that are not subexpressions match
580                 if (!match_same_type(ex_to<basic>(pattern)))
581                         return false;
582
583                 // Otherwise the subexpressions must match one-to-one
584                 for (size_t i=0; i<nops(); i++)
585                         if (!op(i).match(pattern.op(i), repl_lst))
586                                 return false;
587
588                 // Looks similar enough, match found
589                 return true;
590         }
591 }
592
593 /** Helper function for subs(). Does not recurse into subexpressions. */
594 ex basic::subs_one_level(const exmap & m, unsigned options) const
595 {
596         exmap::const_iterator it;
597
598         if (options & subs_options::no_pattern) {
599                 it = m.find(*this);
600                 if (it != m.end())
601                         return it->second;
602         } else {
603                 for (it = m.begin(); it != m.end(); ++it) {
604                         lst repl_lst;
605                         if (match(ex_to<basic>(it->first), repl_lst))
606                                 return it->second.subs(repl_lst, options | subs_options::no_pattern); // avoid infinite recursion when re-substituting the wildcards
607                 }
608         }
609
610         return *this;
611 }
612
613 /** Substitute a set of objects by arbitrary expressions. The ex returned
614  *  will already be evaluated. */
615 ex basic::subs(const exmap & m, unsigned options) const
616 {
617         size_t num = nops();
618         if (num) {
619
620                 // Substitute in subexpressions
621                 for (size_t i=0; i<num; i++) {
622                         const ex & orig_op = op(i);
623                         const ex & subsed_op = orig_op.subs(m, options);
624                         if (!are_ex_trivially_equal(orig_op, subsed_op)) {
625
626                                 // Something changed, clone the object
627                                 basic *copy = duplicate();
628                                 copy->setflag(status_flags::dynallocated);
629                                 copy->clearflag(status_flags::hash_calculated | status_flags::expanded);
630
631                                 // Substitute the changed operand
632                                 copy->let_op(i++) = subsed_op;
633
634                                 // Substitute the other operands
635                                 for (; i<num; i++)
636                                         copy->let_op(i) = op(i).subs(m, options);
637
638                                 // Perform substitutions on the new object as a whole
639                                 return copy->subs_one_level(m, options);
640                         }
641                 }
642         }
643
644         // Nothing changed or no subexpressions
645         return subs_one_level(m, options);
646 }
647
648 /** Default interface of nth derivative ex::diff(s, n).  It should be called
649  *  instead of ::derivative(s) for first derivatives and for nth derivatives it
650  *  just recurses down.
651  *
652  *  @param s symbol to differentiate in
653  *  @param nth order of differentiation
654  *  @see ex::diff */
655 ex basic::diff(const symbol & s, unsigned nth) const
656 {
657         // trivial: zeroth derivative
658         if (nth==0)
659                 return ex(*this);
660         
661         // evaluate unevaluated *this before differentiating
662         if (!(flags & status_flags::evaluated))
663                 return ex(*this).diff(s, nth);
664         
665         ex ndiff = this->derivative(s);
666         while (!ndiff.is_zero() &&    // stop differentiating zeros
667                nth>1) {
668                 ndiff = ndiff.diff(s);
669                 --nth;
670         }
671         return ndiff;
672 }
673
674 /** Return a vector containing the free indices of an expression. */
675 exvector basic::get_free_indices() const
676 {
677         return exvector(); // return an empty exvector
678 }
679
680 ex basic::eval_ncmul(const exvector & v) const
681 {
682         return hold_ncmul(v);
683 }
684
685 // protected
686
687 /** Function object to be applied by basic::derivative(). */
688 struct derivative_map_function : public map_function {
689         const symbol &s;
690         derivative_map_function(const symbol &sym) : s(sym) {}
691         ex operator()(const ex & e) { return diff(e, s); }
692 };
693
694 /** Default implementation of ex::diff(). It maps the operation on the
695  *  operands (or returns 0 when the object has no operands).
696  *
697  *  @see ex::diff */
698 ex basic::derivative(const symbol & s) const
699 {
700         if (nops() == 0)
701                 return _ex0;
702         else {
703                 derivative_map_function map_derivative(s);
704                 return map(map_derivative);
705         }
706 }
707
708 /** Returns order relation between two objects of same type.  This needs to be
709  *  implemented by each class. It may never return anything else than 0,
710  *  signalling equality, or +1 and -1 signalling inequality and determining
711  *  the canonical ordering.  (Perl hackers will wonder why C++ doesn't feature
712  *  the spaceship operator <=> for denoting just this.) */
713 int basic::compare_same_type(const basic & other) const
714 {
715         return compare_pointers(this, &other);
716 }
717
718 /** Returns true if two objects of same type are equal.  Normally needs
719  *  not be reimplemented as long as it wasn't overwritten by some parent
720  *  class, since it just calls compare_same_type().  The reason why this
721  *  function exists is that sometimes it is easier to determine equality
722  *  than an order relation and then it can be overridden. */
723 bool basic::is_equal_same_type(const basic & other) const
724 {
725         return compare_same_type(other)==0;
726 }
727
728 /** Returns true if the attributes of two objects are similar enough for
729  *  a match. This function must not match subexpressions (this is already
730  *  done by basic::match()). Only attributes not accessible by op() should
731  *  be compared. This is also the reason why this function doesn't take the
732  *  wildcard replacement list from match() as an argument: only subexpressions
733  *  are subject to wildcard matches. Also, this function only needs to be
734  *  implemented for container classes because is_equal_same_type() is
735  *  automatically used instead of match_same_type() if nops() == 0.
736  *
737  *  @see basic::match */
738 bool basic::match_same_type(const basic & other) const
739 {
740         // The default is to only consider subexpressions, but not any other
741         // attributes
742         return true;
743 }
744
745 unsigned basic::return_type() const
746 {
747         return return_types::commutative;
748 }
749
750 unsigned basic::return_type_tinfo() const
751 {
752         return tinfo();
753 }
754
755 /** Compute the hash value of an object and if it makes sense to store it in
756  *  the objects status_flags, do so.  The method inherited from class basic
757  *  computes a hash value based on the type and hash values of possible
758  *  members.  For this reason it is well suited for container classes but
759  *  atomic classes should override this implementation because otherwise they
760  *  would all end up with the same hashvalue. */
761 unsigned basic::calchash() const
762 {
763         unsigned v = golden_ratio_hash(tinfo());
764         for (size_t i=0; i<nops(); i++) {
765                 v = rotate_left(v);
766                 v ^= this->op(i).gethash();
767         }
768
769         // store calculated hash value only if object is already evaluated
770         if (flags & status_flags::evaluated) {
771                 setflag(status_flags::hash_calculated);
772                 hashvalue = v;
773         }
774
775         return v;
776 }
777
778 /** Function object to be applied by basic::expand(). */
779 struct expand_map_function : public map_function {
780         unsigned options;
781         expand_map_function(unsigned o) : options(o) {}
782         ex operator()(const ex & e) { return expand(e, options); }
783 };
784
785 /** Expand expression, i.e. multiply it out and return the result as a new
786  *  expression. */
787 ex basic::expand(unsigned options) const
788 {
789         if (nops() == 0)
790                 return (options == 0) ? setflag(status_flags::expanded) : *this;
791         else {
792                 expand_map_function map_expand(options);
793                 return ex_to<basic>(map(map_expand)).setflag(options == 0 ? status_flags::expanded : 0);
794         }
795 }
796
797
798 //////////
799 // non-virtual functions in this class
800 //////////
801
802 // public
803
804 /** Compare objects syntactically to establish canonical ordering.
805  *  All compare functions return: -1 for *this less than other, 0 equal,
806  *  1 greater. */
807 int basic::compare(const basic & other) const
808 {
809         const unsigned hash_this = gethash();
810         const unsigned hash_other = other.gethash();
811         if (hash_this<hash_other) return -1;
812         if (hash_this>hash_other) return 1;
813
814         const unsigned typeid_this = tinfo();
815         const unsigned typeid_other = other.tinfo();
816         if (typeid_this==typeid_other) {
817                 GINAC_ASSERT(typeid(*this)==typeid(other));
818 //              int cmpval = compare_same_type(other);
819 //              if (cmpval!=0) {
820 //                      std::cout << "hash collision, same type: " 
821 //                                << *this << " and " << other << std::endl;
822 //                      this->print(print_tree(std::cout));
823 //                      std::cout << " and ";
824 //                      other.print(print_tree(std::cout));
825 //                      std::cout << std::endl;
826 //              }
827 //              return cmpval;
828                 return compare_same_type(other);
829         } else {
830 //              std::cout << "hash collision, different types: " 
831 //                        << *this << " and " << other << std::endl;
832 //              this->print(print_tree(std::cout));
833 //              std::cout << " and ";
834 //              other.print(print_tree(std::cout));
835 //              std::cout << std::endl;
836                 return (typeid_this<typeid_other ? -1 : 1);
837         }
838 }
839
840 /** Test for syntactic equality.
841  *  This is only a quick test, meaning objects should be in the same domain.
842  *  You might have to .expand(), .normal() objects first, depending on the
843  *  domain of your computation, to get a more reliable answer.
844  *
845  *  @see is_equal_same_type */
846 bool basic::is_equal(const basic & other) const
847 {
848         if (this->gethash()!=other.gethash())
849                 return false;
850         if (this->tinfo()!=other.tinfo())
851                 return false;
852         
853         GINAC_ASSERT(typeid(*this)==typeid(other));
854         
855         return is_equal_same_type(other);
856 }
857
858 // protected
859
860 /** Stop further evaluation.
861  *
862  *  @see basic::eval */
863 const basic & basic::hold() const
864 {
865         return setflag(status_flags::evaluated);
866 }
867
868 /** Ensure the object may be modified without hurting others, throws if this
869  *  is not the case. */
870 void basic::ensure_if_modifiable() const
871 {
872         if (refcount > 1)
873                 throw(std::runtime_error("cannot modify multiply referenced object"));
874         clearflag(status_flags::hash_calculated | status_flags::evaluated);
875 }
876
877 //////////
878 // global variables
879 //////////
880
881 int max_recursion_level = 1024;
882
883 } // namespace GiNaC