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