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