3 * Implementation of GiNaC's ABC. */
6 * GiNaC Copyright (C) 1999-2008 Johannes Gutenberg University Mainz, Germany
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.
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.
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
35 #include "relational.h"
36 #include "operators.h"
44 GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(basic, void,
45 print_func<print_context>(&basic::do_print).
46 print_func<print_tree>(&basic::do_print_tree).
47 print_func<print_python_repr>(&basic::do_print_python_repr))
50 // default constructor, destructor, copy constructor and assignment operator
55 /** basic copy constructor: implicitly assumes that the other class is of
56 * the exact same type (as it's used by duplicate()), so it can copy the
57 * tinfo_key and the hash value. */
58 basic::basic(const basic & other) : tinfo_key(other.tinfo_key), flags(other.flags & ~status_flags::dynallocated), hashvalue(other.hashvalue)
62 /** basic assignment operator: the other object might be of a derived class. */
63 const basic & basic::operator=(const basic & other)
65 unsigned fl = other.flags & ~status_flags::dynallocated;
66 if (typeid(*this) != typeid(other)) {
67 // The other object is of a derived class, so clear the flags as they
68 // might no longer apply (especially hash_calculated). Oh, and don't
69 // copy the tinfo_key: it is already set correctly for this object.
70 fl &= ~(status_flags::evaluated | status_flags::expanded | status_flags::hash_calculated);
72 // The objects are of the exact same class, so copy the hash value.
73 hashvalue = other.hashvalue;
94 /** Construct object from archive_node. */
95 basic::basic(const archive_node &n, lst &sym_lst) : flags(0)
97 // Reconstruct tinfo_key from class name
98 std::string class_name;
99 if (n.find_string("class", class_name))
100 tinfo_key = find_tinfo_key(class_name);
102 throw (std::runtime_error("archive node contains no class name"));
105 /** Unarchive the object. */
106 DEFAULT_UNARCHIVE(basic)
108 /** Archive the object. */
109 void basic::archive(archive_node &n) const
111 n.add_string("class", class_name());
115 // new virtual functions which can be overridden by derived classes
120 /** Output to stream. This performs double dispatch on the dynamic type of
121 * *this and the dynamic type of the supplied print context.
122 * @param c print context object that describes the output formatting
123 * @param level value that is used to identify the precedence or indentation
124 * level for placing parentheses and formatting */
125 void basic::print(const print_context & c, unsigned level) const
127 print_dispatch(get_class_info(), c, level);
130 /** Like print(), but dispatch to the specified class. Can be used by
131 * implementations of print methods to dispatch to the method of the
134 * @see basic::print */
135 void basic::print_dispatch(const registered_class_info & ri, const print_context & c, unsigned level) const
137 // Double dispatch on object type and print_context type
138 const registered_class_info * reg_info = &ri;
139 const print_context_class_info * pc_info = &c.get_class_info();
142 const std::vector<print_functor> & pdt = reg_info->options.get_print_dispatch_table();
145 unsigned id = pc_info->options.get_id();
146 if (id >= pdt.size() || !(pdt[id].is_valid())) {
148 // Method not found, try parent print_context class
149 const print_context_class_info * parent_pc_info = pc_info->get_parent();
150 if (parent_pc_info) {
151 pc_info = parent_pc_info;
155 // Method still not found, try parent class
156 const registered_class_info * parent_reg_info = reg_info->get_parent();
157 if (parent_reg_info) {
158 reg_info = parent_reg_info;
159 pc_info = &c.get_class_info();
163 // Method still not found. This shouldn't happen because basic (the
164 // base class of the algebraic hierarchy) registers a method for
165 // print_context (the base class of the print context hierarchy),
166 // so if we end up here, there's something wrong with the class
168 throw (std::runtime_error(std::string("basic::print(): method for ") + class_name() + "/" + c.class_name() + " not found"));
173 pdt[id](*this, c, level);
177 /** Default output to stream. */
178 void basic::do_print(const print_context & c, unsigned level) const
180 c.s << "[" << class_name() << " object]";
183 /** Tree output to stream. */
184 void basic::do_print_tree(const print_tree & c, unsigned level) const
186 c.s << std::string(level, ' ') << class_name() << " @" << this
187 << std::hex << ", hash=0x" << hashvalue << ", flags=0x" << flags << std::dec;
189 c.s << ", nops=" << nops();
191 for (size_t i=0; i<nops(); ++i)
192 op(i).print(c, level + c.delta_indent);
195 /** Python parsable output to stream. */
196 void basic::do_print_python_repr(const print_python_repr & c, unsigned level) const
198 c.s << class_name() << "()";
201 /** Little wrapper around print to be called within a debugger.
202 * This is needed because you cannot call foo.print(cout) from within the
203 * debugger because it might not know what cout is. This method can be
204 * invoked with no argument and it will simply print to stdout.
207 * @see basic::dbgprinttree */
208 void basic::dbgprint() const
210 this->print(print_dflt(std::cerr));
211 std::cerr << std::endl;
214 /** Little wrapper around printtree to be called within a debugger.
216 * @see basic::dbgprint */
217 void basic::dbgprinttree() const
219 this->print(print_tree(std::cerr));
222 /** Return relative operator precedence (for parenthezing output). */
223 unsigned basic::precedence() const
228 /** Information about the object.
230 * @see class info_flags */
231 bool basic::info(unsigned inf) const
233 // all possible properties are false for basic objects
237 /** Number of operands/members. */
238 size_t basic::nops() const
240 // iterating from 0 to nops() on atomic objects should be an empty loop,
241 // and accessing their elements is a range error. Container objects should
246 /** Return operand/member at position i. */
247 ex basic::op(size_t i) const
249 throw(std::range_error(std::string("basic::op(): ") + class_name() + std::string(" has no operands")));
252 /** Return modifyable operand/member at position i. */
253 ex & basic::let_op(size_t i)
255 ensure_if_modifiable();
256 throw(std::range_error(std::string("basic::let_op(): ") + class_name() + std::string(" has no operands")));
259 ex basic::operator[](const ex & index) const
261 if (is_exactly_a<numeric>(index))
262 return op(static_cast<size_t>(ex_to<numeric>(index).to_int()));
264 throw(std::invalid_argument(std::string("non-numeric indices not supported by ") + class_name()));
267 ex basic::operator[](size_t i) const
272 ex & basic::operator[](const ex & index)
274 if (is_exactly_a<numeric>(index))
275 return let_op(ex_to<numeric>(index).to_int());
277 throw(std::invalid_argument(std::string("non-numeric indices not supported by ") + class_name()));
280 ex & basic::operator[](size_t i)
285 /** Test for occurrence of a pattern. An object 'has' a pattern if it matches
286 * the pattern itself or one of the children 'has' it. As a consequence
287 * (according to the definition of children) given e=x+y+z, e.has(x) is true
288 * but e.has(x+y) is false. */
289 bool basic::has(const ex & pattern, unsigned options) const
292 if (match(pattern, repl_lst))
294 for (size_t i=0; i<nops(); i++)
295 if (op(i).has(pattern, options))
301 /** Construct new expression by applying the specified function to all
302 * sub-expressions (one level only, not recursively). */
303 ex basic::map(map_function & f) const
310 for (size_t i=0; i<num; i++) {
311 const ex & o = op(i);
313 if (!are_ex_trivially_equal(o, n)) {
321 copy->setflag(status_flags::dynallocated);
322 copy->clearflag(status_flags::hash_calculated | status_flags::expanded);
328 /** Check whether this is a polynomial in the given variables. */
329 bool basic::is_polynomial(const ex & var) const
331 return !has(var) || is_equal(ex_to<basic>(var));
334 /** Return degree of highest power in object s. */
335 int basic::degree(const ex & s) const
337 return is_equal(ex_to<basic>(s)) ? 1 : 0;
340 /** Return degree of lowest power in object s. */
341 int basic::ldegree(const ex & s) const
343 return is_equal(ex_to<basic>(s)) ? 1 : 0;
346 /** Return coefficient of degree n in object s. */
347 ex basic::coeff(const ex & s, int n) const
349 if (is_equal(ex_to<basic>(s)))
350 return n==1 ? _ex1 : _ex0;
352 return n==0 ? *this : _ex0;
355 /** Sort expanded expression in terms of powers of some object(s).
356 * @param s object(s) to sort in
357 * @param distributed recursive or distributed form (only used when s is a list) */
358 ex basic::collect(const ex & s, bool distributed) const
363 // List of objects specified
367 return collect(s.op(0));
369 else if (distributed) {
374 const lst& l(ex_to<lst>(s));
378 for (const_iterator xi=x.begin(); xi!=x.end(); ++xi) {
381 for (lst::const_iterator li=l.begin(); li!=l.end(); ++li) {
382 int cexp = pre_coeff.degree(*li);
383 pre_coeff = pre_coeff.coeff(*li, cexp);
384 key *= pow(*li, cexp);
386 exmap::iterator ci = cmap.find(key);
387 if (ci != cmap.end())
388 ci->second += pre_coeff;
390 cmap.insert(exmap::value_type(key, pre_coeff));
394 for (exmap::const_iterator mi=cmap.begin(); mi != cmap.end(); ++mi)
395 resv.push_back((mi->first)*(mi->second));
396 return (new add(resv))->setflag(status_flags::dynallocated);
402 size_t n = s.nops() - 1;
413 // Only one object specified
414 for (int n=this->ldegree(s); n<=this->degree(s); ++n)
415 x += this->coeff(s,n)*power(s,n);
418 // correct for lost fractional arguments and return
419 return x + (*this - x).expand();
422 /** Perform automatic non-interruptive term rewriting rules. */
423 ex basic::eval(int level) const
425 // There is nothing to do for basic objects:
429 /** Function object to be applied by basic::evalf(). */
430 struct evalf_map_function : public map_function {
432 evalf_map_function(int l) : level(l) {}
433 ex operator()(const ex & e) { return evalf(e, level); }
436 /** Evaluate object numerically. */
437 ex basic::evalf(int level) const
444 else if (level == -max_recursion_level)
445 throw(std::runtime_error("max recursion level reached"));
447 evalf_map_function map_evalf(level - 1);
448 return map(map_evalf);
453 /** Function object to be applied by basic::evalm(). */
454 struct evalm_map_function : public map_function {
455 ex operator()(const ex & e) { return evalm(e); }
458 /** Evaluate sums, products and integer powers of matrices. */
459 ex basic::evalm() const
464 return map(map_evalm);
467 /** Function object to be applied by basic::eval_integ(). */
468 struct eval_integ_map_function : public map_function {
469 ex operator()(const ex & e) { return eval_integ(e); }
472 /** Evaluate integrals, if result is known. */
473 ex basic::eval_integ() const
478 return map(map_eval_integ);
481 /** Perform automatic symbolic evaluations on indexed expression that
482 * contains this object as the base expression. */
483 ex basic::eval_indexed(const basic & i) const
484 // this function can't take a "const ex & i" because that would result
485 // in an infinite eval() loop
487 // There is nothing to do for basic objects
491 /** Add two indexed expressions. They are guaranteed to be of class indexed
492 * (or a subclass) and their indices are compatible. This function is used
493 * internally by simplify_indexed().
495 * @param self First indexed expression; its base object is *this
496 * @param other Second indexed expression
497 * @return sum of self and other
498 * @see ex::simplify_indexed() */
499 ex basic::add_indexed(const ex & self, const ex & other) const
504 /** Multiply an indexed expression with a scalar. This function is used
505 * internally by simplify_indexed().
507 * @param self Indexed expression; its base object is *this
508 * @param other Numeric value
509 * @return product of self and other
510 * @see ex::simplify_indexed() */
511 ex basic::scalar_mul_indexed(const ex & self, const numeric & other) const
516 /** Try to contract two indexed expressions that appear in the same product.
517 * If a contraction exists, the function overwrites one or both of the
518 * expressions and returns true. Otherwise it returns false. It is
519 * guaranteed that both expressions are of class indexed (or a subclass)
520 * and that at least one dummy index has been found. This functions is
521 * used internally by simplify_indexed().
523 * @param self Pointer to first indexed expression; its base object is *this
524 * @param other Pointer to second indexed expression
525 * @param v The complete vector of factors
526 * @return true if the contraction was successful, false otherwise
527 * @see ex::simplify_indexed() */
528 bool basic::contract_with(exvector::iterator self, exvector::iterator other, exvector & v) const
534 /** Check whether the expression matches a given pattern. For every wildcard
535 * object in the pattern, a pair with the wildcard as a key and matching
536 * expression as a value is added to repl_lst. */
537 bool basic::match(const ex & pattern, exmap& repl_lst) const
540 Sweet sweet shapes, sweet sweet shapes,
541 That's the key thing, right right.
542 Feed feed face, feed feed shapes,
543 But who is the king tonight?
544 Who is the king tonight?
545 Pattern is the thing, the key thing-a-ling,
546 But who is the king of Pattern?
547 But who is the king, the king thing-a-ling,
548 Who is the king of Pattern?
549 Bog is the king, the king thing-a-ling,
550 Bog is the king of Pattern.
551 Ba bu-bu-bu-bu bu-bu-bu-bu-bu-bu bu-bu
552 Bog is the king of Pattern.
555 if (is_exactly_a<wildcard>(pattern)) {
557 // Wildcard matches anything, but check whether we already have found
558 // a match for that wildcard first (if so, the earlier match must be
559 // the same expression)
560 for (exmap::const_iterator it = repl_lst.begin(); it != repl_lst.end(); ++it) {
561 if (it->first.is_equal(pattern))
562 return is_equal(ex_to<basic>(it->second));
564 repl_lst[pattern] = *this;
569 // Expression must be of the same type as the pattern
570 if (typeid(*this) != typeid(ex_to<basic>(pattern)))
573 // Number of subexpressions must match
574 if (nops() != pattern.nops())
577 // No subexpressions? Then just compare the objects (there can't be
578 // wildcards in the pattern)
580 return is_equal_same_type(ex_to<basic>(pattern));
582 // Check whether attributes that are not subexpressions match
583 if (!match_same_type(ex_to<basic>(pattern)))
586 // Even if the expression does not match the pattern, some of
587 // its subexpressions could match it. For example, x^5*y^(-1)
588 // does not match the pattern $0^5, but its subexpression x^5
589 // does. So, save repl_lst in order to not add bogus entries.
590 exmap tmp_repl = repl_lst;
591 // Otherwise the subexpressions must match one-to-one
592 for (size_t i=0; i<nops(); i++)
593 if (!op(i).match(pattern.op(i), tmp_repl))
596 // Looks similar enough, match found
602 /** Helper function for subs(). Does not recurse into subexpressions. */
603 ex basic::subs_one_level(const exmap & m, unsigned options) const
605 exmap::const_iterator it;
607 if (options & subs_options::no_pattern) {
614 for (it = m.begin(); it != m.end(); ++it) {
616 if (match(ex_to<basic>(it->first), repl_lst))
617 return it->second.subs(repl_lst, options | subs_options::no_pattern);
618 // avoid infinite recursion when re-substituting the wildcards
625 /** Substitute a set of objects by arbitrary expressions. The ex returned
626 * will already be evaluated. */
627 ex basic::subs(const exmap & m, unsigned options) const
632 // Substitute in subexpressions
633 for (size_t i=0; i<num; i++) {
634 const ex & orig_op = op(i);
635 const ex & subsed_op = orig_op.subs(m, options);
636 if (!are_ex_trivially_equal(orig_op, subsed_op)) {
638 // Something changed, clone the object
639 basic *copy = duplicate();
640 copy->setflag(status_flags::dynallocated);
641 copy->clearflag(status_flags::hash_calculated | status_flags::expanded);
643 // Substitute the changed operand
644 copy->let_op(i++) = subsed_op;
646 // Substitute the other operands
648 copy->let_op(i) = op(i).subs(m, options);
650 // Perform substitutions on the new object as a whole
651 return copy->subs_one_level(m, options);
656 // Nothing changed or no subexpressions
657 return subs_one_level(m, options);
660 /** Default interface of nth derivative ex::diff(s, n). It should be called
661 * instead of ::derivative(s) for first derivatives and for nth derivatives it
662 * just recurses down.
664 * @param s symbol to differentiate in
665 * @param nth order of differentiation
667 ex basic::diff(const symbol & s, unsigned nth) const
669 // trivial: zeroth derivative
673 // evaluate unevaluated *this before differentiating
674 if (!(flags & status_flags::evaluated))
675 return ex(*this).diff(s, nth);
677 ex ndiff = this->derivative(s);
678 while (!ndiff.is_zero() && // stop differentiating zeros
680 ndiff = ndiff.diff(s);
686 /** Return a vector containing the free indices of an expression. */
687 exvector basic::get_free_indices() const
689 return exvector(); // return an empty exvector
692 ex basic::conjugate() const
697 ex basic::real_part() const
699 return real_part_function(*this).hold();
702 ex basic::imag_part() const
704 return imag_part_function(*this).hold();
707 ex basic::eval_ncmul(const exvector & v) const
709 return hold_ncmul(v);
714 /** Function object to be applied by basic::derivative(). */
715 struct derivative_map_function : public map_function {
717 derivative_map_function(const symbol &sym) : s(sym) {}
718 ex operator()(const ex & e) { return diff(e, s); }
721 /** Default implementation of ex::diff(). It maps the operation on the
722 * operands (or returns 0 when the object has no operands).
725 ex basic::derivative(const symbol & s) const
730 derivative_map_function map_derivative(s);
731 return map(map_derivative);
735 /** Returns order relation between two objects of same type. This needs to be
736 * implemented by each class. It may never return anything else than 0,
737 * signalling equality, or +1 and -1 signalling inequality and determining
738 * the canonical ordering. (Perl hackers will wonder why C++ doesn't feature
739 * the spaceship operator <=> for denoting just this.) */
740 int basic::compare_same_type(const basic & other) const
742 return compare_pointers(this, &other);
745 /** Returns true if two objects of same type are equal. Normally needs
746 * not be reimplemented as long as it wasn't overwritten by some parent
747 * class, since it just calls compare_same_type(). The reason why this
748 * function exists is that sometimes it is easier to determine equality
749 * than an order relation and then it can be overridden. */
750 bool basic::is_equal_same_type(const basic & other) const
752 return compare_same_type(other)==0;
755 /** Returns true if the attributes of two objects are similar enough for
756 * a match. This function must not match subexpressions (this is already
757 * done by basic::match()). Only attributes not accessible by op() should
758 * be compared. This is also the reason why this function doesn't take the
759 * wildcard replacement list from match() as an argument: only subexpressions
760 * are subject to wildcard matches. Also, this function only needs to be
761 * implemented for container classes because is_equal_same_type() is
762 * automatically used instead of match_same_type() if nops() == 0.
764 * @see basic::match */
765 bool basic::match_same_type(const basic & other) const
767 // The default is to only consider subexpressions, but not any other
772 unsigned basic::return_type() const
774 return return_types::commutative;
777 tinfo_t basic::return_type_tinfo() const
782 /** Compute the hash value of an object and if it makes sense to store it in
783 * the objects status_flags, do so. The method inherited from class basic
784 * computes a hash value based on the type and hash values of possible
785 * members. For this reason it is well suited for container classes but
786 * atomic classes should override this implementation because otherwise they
787 * would all end up with the same hashvalue. */
788 unsigned basic::calchash() const
790 const void* this_tinfo = (const void*)typeid(*this).name();
791 unsigned v = golden_ratio_hash((p_int)this_tinfo);
792 for (size_t i=0; i<nops(); i++) {
794 v ^= this->op(i).gethash();
797 // store calculated hash value only if object is already evaluated
798 if (flags & status_flags::evaluated) {
799 setflag(status_flags::hash_calculated);
806 /** Function object to be applied by basic::expand(). */
807 struct expand_map_function : public map_function {
809 expand_map_function(unsigned o) : options(o) {}
810 ex operator()(const ex & e) { return e.expand(options); }
813 /** Expand expression, i.e. multiply it out and return the result as a new
815 ex basic::expand(unsigned options) const
818 return (options == 0) ? setflag(status_flags::expanded) : *this;
820 expand_map_function map_expand(options);
821 return ex_to<basic>(map(map_expand)).setflag(options == 0 ? status_flags::expanded : 0);
827 // non-virtual functions in this class
832 /** Compare objects syntactically to establish canonical ordering.
833 * All compare functions return: -1 for *this less than other, 0 equal,
835 int basic::compare(const basic & other) const
837 #ifdef GINAC_COMPARE_STATISTICS
838 compare_statistics.total_basic_compares++;
840 const unsigned hash_this = gethash();
841 const unsigned hash_other = other.gethash();
842 if (hash_this<hash_other) return -1;
843 if (hash_this>hash_other) return 1;
844 #ifdef GINAC_COMPARE_STATISTICS
845 compare_statistics.compare_same_hashvalue++;
848 const std::type_info& typeid_this = typeid(*this);
849 const std::type_info& typeid_other = typeid(other);
850 if (typeid_this == typeid_other) {
851 // int cmpval = compare_same_type(other);
853 // std::cout << "hash collision, same type: "
854 // << *this << " and " << other << std::endl;
855 // this->print(print_tree(std::cout));
856 // std::cout << " and ";
857 // other.print(print_tree(std::cout));
858 // std::cout << std::endl;
861 #ifdef GINAC_COMPARE_STATISTICS
862 compare_statistics.compare_same_type++;
864 return compare_same_type(other);
866 // std::cout << "hash collision, different types: "
867 // << *this << " and " << other << std::endl;
868 // this->print(print_tree(std::cout));
869 // std::cout << " and ";
870 // other.print(print_tree(std::cout));
871 // std::cout << std::endl;
872 return (typeid_this.before(typeid_other) ? -1 : 1);
876 /** Test for syntactic equality.
877 * This is only a quick test, meaning objects should be in the same domain.
878 * You might have to .expand(), .normal() objects first, depending on the
879 * domain of your computation, to get a more reliable answer.
881 * @see is_equal_same_type */
882 bool basic::is_equal(const basic & other) const
884 #ifdef GINAC_COMPARE_STATISTICS
885 compare_statistics.total_basic_is_equals++;
887 if (this->gethash()!=other.gethash())
889 #ifdef GINAC_COMPARE_STATISTICS
890 compare_statistics.is_equal_same_hashvalue++;
892 if (typeid(*this) != typeid(other))
895 #ifdef GINAC_COMPARE_STATISTICS
896 compare_statistics.is_equal_same_type++;
898 return is_equal_same_type(other);
903 /** Stop further evaluation.
905 * @see basic::eval */
906 const basic & basic::hold() const
908 return setflag(status_flags::evaluated);
911 /** Ensure the object may be modified without hurting others, throws if this
912 * is not the case. */
913 void basic::ensure_if_modifiable() const
915 if (get_refcount() > 1)
916 throw(std::runtime_error("cannot modify multiply referenced object"));
917 clearflag(status_flags::hash_calculated | status_flags::evaluated);
924 int max_recursion_level = 1024;
927 #ifdef GINAC_COMPARE_STATISTICS
928 compare_statistics_t::~compare_statistics_t()
930 std::clog << "ex::compare() called " << total_compares << " times" << std::endl;
931 std::clog << "nontrivial compares: " << nontrivial_compares << " times" << std::endl;
932 std::clog << "basic::compare() called " << total_basic_compares << " times" << std::endl;
933 std::clog << "same hashvalue in compare(): " << compare_same_hashvalue << " times" << std::endl;
934 std::clog << "compare_same_type() called " << compare_same_type << " times" << std::endl;
935 std::clog << std::endl;
936 std::clog << "ex::is_equal() called " << total_is_equals << " times" << std::endl;
937 std::clog << "nontrivial is_equals: " << nontrivial_is_equals << " times" << std::endl;
938 std::clog << "basic::is_equal() called " << total_basic_is_equals << " times" << std::endl;
939 std::clog << "same hashvalue in is_equal(): " << is_equal_same_hashvalue << " times" << std::endl;
940 std::clog << "is_equal_same_type() called " << is_equal_same_type << " times" << std::endl;
941 std::clog << std::endl;
942 std::clog << "basic::gethash() called " << total_gethash << " times" << std::endl;
943 std::clog << "used cached hashvalue " << gethash_cached << " times" << std::endl;
946 compare_statistics_t compare_statistics;