- New tinfo mechanism
[ginac.git] / ginac / basic.cpp
1 /** @file basic.cpp
2  *
3  *  Implementation of GiNaC's ABC. */
4
5 /*
6  *  GiNaC Copyright (C) 1999-2006 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 #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)
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 (tinfo_key != other.tinfo_key) {
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) 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 = 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 /** Return degree of highest power in object s. */
329 int basic::degree(const ex & s) const
330 {
331         return is_equal(ex_to<basic>(s)) ? 1 : 0;
332 }
333
334 /** Return degree of lowest power in object s. */
335 int basic::ldegree(const ex & s) const
336 {
337         return is_equal(ex_to<basic>(s)) ? 1 : 0;
338 }
339
340 /** Return coefficient of degree n in object s. */
341 ex basic::coeff(const ex & s, int n) const
342 {
343         if (is_equal(ex_to<basic>(s)))
344                 return n==1 ? _ex1 : _ex0;
345         else
346                 return n==0 ? *this : _ex0;
347 }
348
349 /** Sort expanded expression in terms of powers of some object(s).
350  *  @param s object(s) to sort in
351  *  @param distributed recursive or distributed form (only used when s is a list) */
352 ex basic::collect(const ex & s, bool distributed) const
353 {
354         ex x;
355         if (is_a<lst>(s)) {
356
357                 // List of objects specified
358                 if (s.nops() == 0)
359                         return *this;
360                 if (s.nops() == 1)
361                         return collect(s.op(0));
362
363                 else if (distributed) {
364
365                         // Get lower/upper degree of all symbols in list
366                         size_t num = s.nops();
367                         struct sym_info {
368                                 ex sym;
369                                 int ldeg, deg;
370                                 int cnt;  // current degree, 'counter'
371                                 ex coeff; // coefficient for degree 'cnt'
372                         };
373                         sym_info *si = new sym_info[num];
374                         ex c = *this;
375                         for (size_t i=0; i<num; i++) {
376                                 si[i].sym = s.op(i);
377                                 si[i].ldeg = si[i].cnt = this->ldegree(si[i].sym);
378                                 si[i].deg = this->degree(si[i].sym);
379                                 c = si[i].coeff = c.coeff(si[i].sym, si[i].cnt);
380                         }
381
382                         while (true) {
383
384                                 // Calculate coeff*x1^c1*...*xn^cn
385                                 ex y = _ex1;
386                                 for (size_t i=0; i<num; i++) {
387                                         int cnt = si[i].cnt;
388                                         y *= power(si[i].sym, cnt);
389                                 }
390                                 x += y * si[num - 1].coeff;
391
392                                 // Increment counters
393                                 size_t n = num - 1;
394                                 while (true) {
395                                         ++si[n].cnt;
396                                         if (si[n].cnt <= si[n].deg) {
397                                                 // Update coefficients
398                                                 ex c;
399                                                 if (n == 0)
400                                                         c = *this;
401                                                 else
402                                                         c = si[n - 1].coeff;
403                                                 for (size_t i=n; i<num; i++)
404                                                         c = si[i].coeff = c.coeff(si[i].sym, si[i].cnt);
405                                                 break;
406                                         }
407                                         if (n == 0)
408                                                 goto done;
409                                         si[n].cnt = si[n].ldeg;
410                                         n--;
411                                 }
412                         }
413
414 done:           delete[] si;
415
416                 } else {
417
418                         // Recursive form
419                         x = *this;
420                         size_t n = s.nops() - 1;
421                         while (true) {
422                                 x = x.collect(s[n]);
423                                 if (n == 0)
424                                         break;
425                                 n--;
426                         }
427                 }
428
429         } else {
430
431                 // Only one object specified
432                 for (int n=this->ldegree(s); n<=this->degree(s); ++n)
433                         x += this->coeff(s,n)*power(s,n);
434         }
435         
436         // correct for lost fractional arguments and return
437         return x + (*this - x).expand();
438 }
439
440 /** Perform automatic non-interruptive term rewriting rules. */
441 ex basic::eval(int level) const
442 {
443         // There is nothing to do for basic objects:
444         return hold();
445 }
446
447 /** Function object to be applied by basic::evalf(). */
448 struct evalf_map_function : public map_function {
449         int level;
450         evalf_map_function(int l) : level(l) {}
451         ex operator()(const ex & e) { return evalf(e, level); }
452 };
453
454 /** Evaluate object numerically. */
455 ex basic::evalf(int level) const
456 {
457         if (nops() == 0)
458                 return *this;
459         else {
460                 if (level == 1)
461                         return *this;
462                 else if (level == -max_recursion_level)
463                         throw(std::runtime_error("max recursion level reached"));
464                 else {
465                         evalf_map_function map_evalf(level - 1);
466                         return map(map_evalf);
467                 }
468         }
469 }
470
471 /** Function object to be applied by basic::evalm(). */
472 struct evalm_map_function : public map_function {
473         ex operator()(const ex & e) { return evalm(e); }
474 } map_evalm;
475
476 /** Evaluate sums, products and integer powers of matrices. */
477 ex basic::evalm() const
478 {
479         if (nops() == 0)
480                 return *this;
481         else
482                 return map(map_evalm);
483 }
484
485 /** Function object to be applied by basic::eval_integ(). */
486 struct eval_integ_map_function : public map_function {
487         ex operator()(const ex & e) { return eval_integ(e); }
488 } map_eval_integ;
489
490 /** Evaluate integrals, if result is known. */
491 ex basic::eval_integ() const
492 {
493         if (nops() == 0)
494                 return *this;
495         else
496                 return map(map_eval_integ);
497 }
498
499 /** Perform automatic symbolic evaluations on indexed expression that
500  *  contains this object as the base expression. */
501 ex basic::eval_indexed(const basic & i) const
502  // this function can't take a "const ex & i" because that would result
503  // in an infinite eval() loop
504 {
505         // There is nothing to do for basic objects
506         return i.hold();
507 }
508
509 /** Add two indexed expressions. They are guaranteed to be of class indexed
510  *  (or a subclass) and their indices are compatible. This function is used
511  *  internally by simplify_indexed().
512  *
513  *  @param self First indexed expression; its base object is *this
514  *  @param other Second indexed expression
515  *  @return sum of self and other 
516  *  @see ex::simplify_indexed() */
517 ex basic::add_indexed(const ex & self, const ex & other) const
518 {
519         return self + other;
520 }
521
522 /** Multiply an indexed expression with a scalar. This function is used
523  *  internally by simplify_indexed().
524  *
525  *  @param self Indexed expression; its base object is *this
526  *  @param other Numeric value
527  *  @return product of self and other
528  *  @see ex::simplify_indexed() */
529 ex basic::scalar_mul_indexed(const ex & self, const numeric & other) const
530 {
531         return self * other;
532 }
533
534 /** Try to contract two indexed expressions that appear in the same product. 
535  *  If a contraction exists, the function overwrites one or both of the
536  *  expressions and returns true. Otherwise it returns false. It is
537  *  guaranteed that both expressions are of class indexed (or a subclass)
538  *  and that at least one dummy index has been found. This functions is
539  *  used internally by simplify_indexed().
540  *
541  *  @param self Pointer to first indexed expression; its base object is *this
542  *  @param other Pointer to second indexed expression
543  *  @param v The complete vector of factors
544  *  @return true if the contraction was successful, false otherwise
545  *  @see ex::simplify_indexed() */
546 bool basic::contract_with(exvector::iterator self, exvector::iterator other, exvector & v) const
547 {
548         // Do nothing
549         return false;
550 }
551
552 /** Check whether the expression matches a given pattern. For every wildcard
553  *  object in the pattern, an expression of the form "wildcard == matching_expression"
554  *  is added to repl_lst. */
555 bool basic::match(const ex & pattern, lst & repl_lst) const
556 {
557 /*
558         Sweet sweet shapes, sweet sweet shapes,
559         That's the key thing, right right.
560         Feed feed face, feed feed shapes,
561         But who is the king tonight?
562         Who is the king tonight?
563         Pattern is the thing, the key thing-a-ling,
564         But who is the king of Pattern?
565         But who is the king, the king thing-a-ling,
566         Who is the king of Pattern?
567         Bog is the king, the king thing-a-ling,
568         Bog is the king of Pattern.
569         Ba bu-bu-bu-bu bu-bu-bu-bu-bu-bu bu-bu
570         Bog is the king of Pattern.
571 */
572
573         if (is_exactly_a<wildcard>(pattern)) {
574
575                 // Wildcard matches anything, but check whether we already have found
576                 // a match for that wildcard first (if so, the earlier match must be
577                 // the same expression)
578                 for (lst::const_iterator it = repl_lst.begin(); it != repl_lst.end(); ++it) {
579                         if (it->op(0).is_equal(pattern))
580                                 return is_equal(ex_to<basic>(it->op(1)));
581                 }
582                 repl_lst.append(pattern == *this);
583                 return true;
584
585         } else {
586
587                 // Expression must be of the same type as the pattern
588                 if (tinfo() != ex_to<basic>(pattern).tinfo())
589                         return false;
590
591                 // Number of subexpressions must match
592                 if (nops() != pattern.nops())
593                         return false;
594
595                 // No subexpressions? Then just compare the objects (there can't be
596                 // wildcards in the pattern)
597                 if (nops() == 0)
598                         return is_equal_same_type(ex_to<basic>(pattern));
599
600                 // Check whether attributes that are not subexpressions match
601                 if (!match_same_type(ex_to<basic>(pattern)))
602                         return false;
603
604                 // Otherwise the subexpressions must match one-to-one
605                 for (size_t i=0; i<nops(); i++)
606                         if (!op(i).match(pattern.op(i), repl_lst))
607                                 return false;
608
609                 // Looks similar enough, match found
610                 return true;
611         }
612 }
613
614 /** Helper function for subs(). Does not recurse into subexpressions. */
615 ex basic::subs_one_level(const exmap & m, unsigned options) const
616 {
617         exmap::const_iterator it;
618
619         if (options & subs_options::no_pattern) {
620                 ex thisex = *this;
621                 it = m.find(thisex);
622                 if (it != m.end())
623                         return it->second;
624                 return thisex;
625         } else {
626                 for (it = m.begin(); it != m.end(); ++it) {
627                         lst repl_lst;
628                         if (match(ex_to<basic>(it->first), repl_lst))
629                                 return it->second.subs(repl_lst, options | subs_options::no_pattern); // avoid infinite recursion when re-substituting the wildcards
630                 }
631         }
632
633         return *this;
634 }
635
636 /** Substitute a set of objects by arbitrary expressions. The ex returned
637  *  will already be evaluated. */
638 ex basic::subs(const exmap & m, unsigned options) const
639 {
640         size_t num = nops();
641         if (num) {
642
643                 // Substitute in subexpressions
644                 for (size_t i=0; i<num; i++) {
645                         const ex & orig_op = op(i);
646                         const ex & subsed_op = orig_op.subs(m, options);
647                         if (!are_ex_trivially_equal(orig_op, subsed_op)) {
648
649                                 // Something changed, clone the object
650                                 basic *copy = duplicate();
651                                 copy->setflag(status_flags::dynallocated);
652                                 copy->clearflag(status_flags::hash_calculated | status_flags::expanded);
653
654                                 // Substitute the changed operand
655                                 copy->let_op(i++) = subsed_op;
656
657                                 // Substitute the other operands
658                                 for (; i<num; i++)
659                                         copy->let_op(i) = op(i).subs(m, options);
660
661                                 // Perform substitutions on the new object as a whole
662                                 return copy->subs_one_level(m, options);
663                         }
664                 }
665         }
666
667         // Nothing changed or no subexpressions
668         return subs_one_level(m, options);
669 }
670
671 /** Default interface of nth derivative ex::diff(s, n).  It should be called
672  *  instead of ::derivative(s) for first derivatives and for nth derivatives it
673  *  just recurses down.
674  *
675  *  @param s symbol to differentiate in
676  *  @param nth order of differentiation
677  *  @see ex::diff */
678 ex basic::diff(const symbol & s, unsigned nth) const
679 {
680         // trivial: zeroth derivative
681         if (nth==0)
682                 return ex(*this);
683         
684         // evaluate unevaluated *this before differentiating
685         if (!(flags & status_flags::evaluated))
686                 return ex(*this).diff(s, nth);
687         
688         ex ndiff = this->derivative(s);
689         while (!ndiff.is_zero() &&    // stop differentiating zeros
690                nth>1) {
691                 ndiff = ndiff.diff(s);
692                 --nth;
693         }
694         return ndiff;
695 }
696
697 /** Return a vector containing the free indices of an expression. */
698 exvector basic::get_free_indices() const
699 {
700         return exvector(); // return an empty exvector
701 }
702
703 ex basic::conjugate() const
704 {
705         return *this;
706 }
707
708 ex basic::eval_ncmul(const exvector & v) const
709 {
710         return hold_ncmul(v);
711 }
712
713 // protected
714
715 /** Function object to be applied by basic::derivative(). */
716 struct derivative_map_function : public map_function {
717         const symbol &s;
718         derivative_map_function(const symbol &sym) : s(sym) {}
719         ex operator()(const ex & e) { return diff(e, s); }
720 };
721
722 /** Default implementation of ex::diff(). It maps the operation on the
723  *  operands (or returns 0 when the object has no operands).
724  *
725  *  @see ex::diff */
726 ex basic::derivative(const symbol & s) const
727 {
728         if (nops() == 0)
729                 return _ex0;
730         else {
731                 derivative_map_function map_derivative(s);
732                 return map(map_derivative);
733         }
734 }
735
736 /** Returns order relation between two objects of same type.  This needs to be
737  *  implemented by each class. It may never return anything else than 0,
738  *  signalling equality, or +1 and -1 signalling inequality and determining
739  *  the canonical ordering.  (Perl hackers will wonder why C++ doesn't feature
740  *  the spaceship operator <=> for denoting just this.) */
741 int basic::compare_same_type(const basic & other) const
742 {
743         return compare_pointers(this, &other);
744 }
745
746 /** Returns true if two objects of same type are equal.  Normally needs
747  *  not be reimplemented as long as it wasn't overwritten by some parent
748  *  class, since it just calls compare_same_type().  The reason why this
749  *  function exists is that sometimes it is easier to determine equality
750  *  than an order relation and then it can be overridden. */
751 bool basic::is_equal_same_type(const basic & other) const
752 {
753         return compare_same_type(other)==0;
754 }
755
756 /** Returns true if the attributes of two objects are similar enough for
757  *  a match. This function must not match subexpressions (this is already
758  *  done by basic::match()). Only attributes not accessible by op() should
759  *  be compared. This is also the reason why this function doesn't take the
760  *  wildcard replacement list from match() as an argument: only subexpressions
761  *  are subject to wildcard matches. Also, this function only needs to be
762  *  implemented for container classes because is_equal_same_type() is
763  *  automatically used instead of match_same_type() if nops() == 0.
764  *
765  *  @see basic::match */
766 bool basic::match_same_type(const basic & other) const
767 {
768         // The default is to only consider subexpressions, but not any other
769         // attributes
770         return true;
771 }
772
773 unsigned basic::return_type() const
774 {
775         return return_types::commutative;
776 }
777
778 const basic* basic::return_type_tinfo() const
779 {
780         return this;
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         unsigned v = golden_ratio_hash((unsigned)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 tinfo_t typeid_this = tinfo();
849         const tinfo_t typeid_other = other.tinfo();
850         if (typeid_this==typeid_other) {
851                 GINAC_ASSERT(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<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 (this->tinfo()!=other.tinfo())
894                 return false;
895         
896         GINAC_ASSERT(typeid(*this)==typeid(other));
897         
898 #ifdef GINAC_COMPARE_STATISTICS
899         compare_statistics.is_equal_same_type++;
900 #endif
901         return is_equal_same_type(other);
902 }
903
904 // protected
905
906 /** Stop further evaluation.
907  *
908  *  @see basic::eval */
909 const basic & basic::hold() const
910 {
911         return setflag(status_flags::evaluated);
912 }
913
914 /** Ensure the object may be modified without hurting others, throws if this
915  *  is not the case. */
916 void basic::ensure_if_modifiable() const
917 {
918         if (get_refcount() > 1)
919                 throw(std::runtime_error("cannot modify multiply referenced object"));
920         clearflag(status_flags::hash_calculated | status_flags::evaluated);
921 }
922
923 //////////
924 // global variables
925 //////////
926
927 int max_recursion_level = 1024;
928
929
930 #ifdef GINAC_COMPARE_STATISTICS
931 compare_statistics_t::~compare_statistics_t()
932 {
933         std::clog << "ex::compare() called " << total_compares << " times" << std::endl;
934         std::clog << "nontrivial compares: " << nontrivial_compares << " times" << std::endl;
935         std::clog << "basic::compare() called " << total_basic_compares << " times" << std::endl;
936         std::clog << "same hashvalue in compare(): " << compare_same_hashvalue << " times" << std::endl;
937         std::clog << "compare_same_type() called " << compare_same_type << " times" << std::endl;
938         std::clog << std::endl;
939         std::clog << "ex::is_equal() called " << total_is_equals << " times" << std::endl;
940         std::clog << "nontrivial is_equals: " << nontrivial_is_equals << " times" << std::endl;
941         std::clog << "basic::is_equal() called " << total_basic_is_equals << " times" << std::endl;
942         std::clog << "same hashvalue in is_equal(): " << is_equal_same_hashvalue << " times" << std::endl;
943         std::clog << "is_equal_same_type() called " << is_equal_same_type << " times" << std::endl;
944         std::clog << std::endl;
945         std::clog << "basic::gethash() called " << total_gethash << " times" << std::endl;
946         std::clog << "used cached hashvalue " << gethash_cached << " times" << std::endl;
947 }
948
949 compare_statistics_t compare_statistics;
950 #endif
951
952 } // namespace GiNaC