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