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