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