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