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