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