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