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