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