595d64500c37435b1ea324b39782d07cc72f88c2
[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 lst & ls, const lst & lr, unsigned options) const
528 {
529         GINAC_ASSERT(ls.nops() == lr.nops());
530
531         lst::const_iterator its, itr;
532
533         if (options & subs_options::subs_no_pattern) {
534                 for (its = ls.begin(), itr = lr.begin(); its != ls.end(); ++its, ++itr) {
535                         if (is_equal(ex_to<basic>(*its)))
536                                 return *itr;
537                 }
538         } else {
539                 for (its = ls.begin(), itr = lr.begin(); its != ls.end(); ++its, ++itr) {
540                         lst repl_lst;
541                         if (match(ex_to<basic>(*its), repl_lst))
542                                 return itr->subs(repl_lst, options | subs_options::subs_no_pattern); // avoid infinite recursion when re-substituting the wildcards
543                 }
544         }
545
546         return *this;
547 }
548
549 /** Substitute a set of objects by arbitrary expressions. The ex returned
550  *  will already be evaluated. */
551 ex basic::subs(const lst & ls, const lst & lr, unsigned options) const
552 {
553         size_t num = nops();
554         if (num) {
555
556                 // Substitute in subexpressions
557                 for (size_t i=0; i<num; i++) {
558                         const ex & orig_op = op(i);
559                         const ex & subsed_op = orig_op.subs(ls, lr, options);
560                         if (!are_ex_trivially_equal(orig_op, subsed_op)) {
561
562                                 // Something changed, clone the object
563                                 basic *copy = duplicate();
564                                 copy->setflag(status_flags::dynallocated);
565                                 copy->clearflag(status_flags::hash_calculated | status_flags::expanded);
566
567                                 // Substitute the changed operand
568                                 copy->let_op(i++) = subsed_op;
569
570                                 // Substitute the other operands
571                                 for (; i<num; i++)
572                                         copy->let_op(i) = op(i).subs(ls, lr, options);
573
574                                 // Perform substitutions on the new object as a whole
575                                 return copy->subs_one_level(ls, lr, options);
576                         }
577                 }
578         }
579
580         // Nothing changed or no subexpressions
581         return subs_one_level(ls, lr, options);
582 }
583
584 /** Default interface of nth derivative ex::diff(s, n).  It should be called
585  *  instead of ::derivative(s) for first derivatives and for nth derivatives it
586  *  just recurses down.
587  *
588  *  @param s symbol to differentiate in
589  *  @param nth order of differentiation
590  *  @see ex::diff */
591 ex basic::diff(const symbol & s, unsigned nth) const
592 {
593         // trivial: zeroth derivative
594         if (nth==0)
595                 return ex(*this);
596         
597         // evaluate unevaluated *this before differentiating
598         if (!(flags & status_flags::evaluated))
599                 return ex(*this).diff(s, nth);
600         
601         ex ndiff = this->derivative(s);
602         while (!ndiff.is_zero() &&    // stop differentiating zeros
603                nth>1) {
604                 ndiff = ndiff.diff(s);
605                 --nth;
606         }
607         return ndiff;
608 }
609
610 /** Return a vector containing the free indices of an expression. */
611 exvector basic::get_free_indices() const
612 {
613         return exvector(); // return an empty exvector
614 }
615
616 ex basic::eval_ncmul(const exvector & v) const
617 {
618         return hold_ncmul(v);
619 }
620
621 // protected
622
623 /** Function object to be applied by basic::derivative(). */
624 struct derivative_map_function : public map_function {
625         const symbol &s;
626         derivative_map_function(const symbol &sym) : s(sym) {}
627         ex operator()(const ex & e) { return diff(e, s); }
628 };
629
630 /** Default implementation of ex::diff(). It maps the operation on the
631  *  operands (or returns 0 when the object has no operands).
632  *
633  *  @see ex::diff */
634 ex basic::derivative(const symbol & s) const
635 {
636         if (nops() == 0)
637                 return _ex0;
638         else {
639                 derivative_map_function map_derivative(s);
640                 return map(map_derivative);
641         }
642 }
643
644 /** Returns order relation between two objects of same type.  This needs to be
645  *  implemented by each class. It may never return anything else than 0,
646  *  signalling equality, or +1 and -1 signalling inequality and determining
647  *  the canonical ordering.  (Perl hackers will wonder why C++ doesn't feature
648  *  the spaceship operator <=> for denoting just this.) */
649 int basic::compare_same_type(const basic & other) const
650 {
651         return compare_pointers(this, &other);
652 }
653
654 /** Returns true if two objects of same type are equal.  Normally needs
655  *  not be reimplemented as long as it wasn't overwritten by some parent
656  *  class, since it just calls compare_same_type().  The reason why this
657  *  function exists is that sometimes it is easier to determine equality
658  *  than an order relation and then it can be overridden. */
659 bool basic::is_equal_same_type(const basic & other) const
660 {
661         return compare_same_type(other)==0;
662 }
663
664 /** Returns true if the attributes of two objects are similar enough for
665  *  a match. This function must not match subexpressions (this is already
666  *  done by basic::match()). Only attributes not accessible by op() should
667  *  be compared. This is also the reason why this function doesn't take the
668  *  wildcard replacement list from match() as an argument: only subexpressions
669  *  are subject to wildcard matches. Also, this function only needs to be
670  *  implemented for container classes because is_equal_same_type() is
671  *  automatically used instead of match_same_type() if nops() == 0.
672  *
673  *  @see basic::match */
674 bool basic::match_same_type(const basic & other) const
675 {
676         // The default is to only consider subexpressions, but not any other
677         // attributes
678         return true;
679 }
680
681 unsigned basic::return_type() const
682 {
683         return return_types::commutative;
684 }
685
686 unsigned basic::return_type_tinfo() const
687 {
688         return tinfo();
689 }
690
691 /** Compute the hash value of an object and if it makes sense to store it in
692  *  the objects status_flags, do so.  The method inherited from class basic
693  *  computes a hash value based on the type and hash values of possible
694  *  members.  For this reason it is well suited for container classes but
695  *  atomic classes should override this implementation because otherwise they
696  *  would all end up with the same hashvalue. */
697 unsigned basic::calchash() const
698 {
699         unsigned v = golden_ratio_hash(tinfo());
700         for (size_t i=0; i<nops(); i++) {
701                 v = rotate_left(v);
702                 v ^= (const_cast<basic *>(this))->op(i).gethash();
703         }
704
705         // store calculated hash value only if object is already evaluated
706         if (flags & status_flags::evaluated) {
707                 setflag(status_flags::hash_calculated);
708                 hashvalue = v;
709         }
710
711         return v;
712 }
713
714 /** Function object to be applied by basic::expand(). */
715 struct expand_map_function : public map_function {
716         unsigned options;
717         expand_map_function(unsigned o) : options(o) {}
718         ex operator()(const ex & e) { return expand(e, options); }
719 };
720
721 /** Expand expression, i.e. multiply it out and return the result as a new
722  *  expression. */
723 ex basic::expand(unsigned options) const
724 {
725         if (nops() == 0)
726                 return (options == 0) ? setflag(status_flags::expanded) : *this;
727         else {
728                 expand_map_function map_expand(options);
729                 return ex_to<basic>(map(map_expand)).setflag(options == 0 ? status_flags::expanded : 0);
730         }
731 }
732
733
734 //////////
735 // non-virtual functions in this class
736 //////////
737
738 // public
739
740 /** Substitute objects in an expression (syntactic substitution) and return
741  *  the result as a new expression.  There are two valid types of
742  *  replacement arguments: 1) a relational like object==ex and 2) a list of
743  *  relationals lst(object1==ex1,object2==ex2,...), which is converted to
744  *  subs(lst(object1,object2,...),lst(ex1,ex2,...)). */
745 ex basic::subs(const ex & e, unsigned options) const
746 {
747         if (e.info(info_flags::relation_equal)) {
748                 return subs(lst(e.lhs()), lst(e.rhs()), options);
749         }
750         if (!e.info(info_flags::list)) {
751                 throw(std::invalid_argument("basic::subs(ex): argument must be a list"));
752         }
753
754         // Split list into two
755         lst ls;
756         lst lr;
757         GINAC_ASSERT(is_a<lst>(e));
758         for (lst::const_iterator it = ex_to<lst>(e).begin(); it != ex_to<lst>(e).end(); ++it) {
759                 ex r = *it;
760                 if (!r.info(info_flags::relation_equal)) {
761                         throw(std::invalid_argument("basic::subs(ex): argument must be a list of equations"));
762                 }
763                 ls.append(r.op(0));
764                 lr.append(r.op(1));
765         }
766         return subs(ls, lr, options);
767 }
768
769 /** Compare objects syntactically to establish canonical ordering.
770  *  All compare functions return: -1 for *this less than other, 0 equal,
771  *  1 greater. */
772 int basic::compare(const basic & other) const
773 {
774         const unsigned hash_this = gethash();
775         const unsigned hash_other = other.gethash();
776         if (hash_this<hash_other) return -1;
777         if (hash_this>hash_other) return 1;
778
779         const unsigned typeid_this = tinfo();
780         const unsigned typeid_other = other.tinfo();
781         if (typeid_this==typeid_other) {
782                 GINAC_ASSERT(typeid(*this)==typeid(other));
783 //              int cmpval = compare_same_type(other);
784 //              if (cmpval!=0) {
785 //                      std::cout << "hash collision, same type: " 
786 //                                << *this << " and " << other << std::endl;
787 //                      this->print(print_tree(std::cout));
788 //                      std::cout << " and ";
789 //                      other.print(print_tree(std::cout));
790 //                      std::cout << std::endl;
791 //              }
792 //              return cmpval;
793                 return compare_same_type(other);
794         } else {
795 //              std::cout << "hash collision, different types: " 
796 //                        << *this << " and " << other << std::endl;
797 //              this->print(print_tree(std::cout));
798 //              std::cout << " and ";
799 //              other.print(print_tree(std::cout));
800 //              std::cout << std::endl;
801                 return (typeid_this<typeid_other ? -1 : 1);
802         }
803 }
804
805 /** Test for syntactic equality.
806  *  This is only a quick test, meaning objects should be in the same domain.
807  *  You might have to .expand(), .normal() objects first, depending on the
808  *  domain of your computation, to get a more reliable answer.
809  *
810  *  @see is_equal_same_type */
811 bool basic::is_equal(const basic & other) const
812 {
813         if (this->gethash()!=other.gethash())
814                 return false;
815         if (this->tinfo()!=other.tinfo())
816                 return false;
817         
818         GINAC_ASSERT(typeid(*this)==typeid(other));
819         
820         return is_equal_same_type(other);
821 }
822
823 // protected
824
825 /** Stop further evaluation.
826  *
827  *  @see basic::eval */
828 const basic & basic::hold() const
829 {
830         return setflag(status_flags::evaluated);
831 }
832
833 /** Ensure the object may be modified without hurting others, throws if this
834  *  is not the case. */
835 void basic::ensure_if_modifiable() const
836 {
837         if (refcount > 1)
838                 throw(std::runtime_error("cannot modify multiply referenced object"));
839         clearflag(status_flags::hash_calculated | status_flags::evaluated);
840 }
841
842 //////////
843 // global variables
844 //////////
845
846 int max_recursion_level = 1024;
847
848 } // namespace GiNaC