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) : 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 throw (std::runtime_error("archive node contains no class name"));
103 /** Unarchive the object. */
104 DEFAULT_UNARCHIVE(basic)
106 /** Archive the object. */
107 void basic::archive(archive_node &n) const
109 n.add_string("class", class_name());
113 // new virtual functions which can be overridden by derived classes
118 /** Output to stream. This performs double dispatch on the dynamic type of
119 * *this and the dynamic type of the supplied print context.
120 * @param c print context object that describes the output formatting
121 * @param level value that is used to identify the precedence or indentation
122 * level for placing parentheses and formatting */
123 void basic::print(const print_context & c, unsigned level) const
125 print_dispatch(get_class_info(), c, level);
128 /** Like print(), but dispatch to the specified class. Can be used by
129 * implementations of print methods to dispatch to the method of the
132 * @see basic::print */
133 void basic::print_dispatch(const registered_class_info & ri, const print_context & c, unsigned level) const
135 // Double dispatch on object type and print_context type
136 const registered_class_info * reg_info = &ri;
137 const print_context_class_info * pc_info = &c.get_class_info();
140 const std::vector<print_functor> & pdt = reg_info->options.get_print_dispatch_table();
143 unsigned id = pc_info->options.get_id();
144 if (id >= pdt.size() || !(pdt[id].is_valid())) {
146 // Method not found, try parent print_context class
147 const print_context_class_info * parent_pc_info = pc_info->get_parent();
148 if (parent_pc_info) {
149 pc_info = parent_pc_info;
153 // Method still not found, try parent class
154 const registered_class_info * parent_reg_info = reg_info->get_parent();
155 if (parent_reg_info) {
156 reg_info = parent_reg_info;
157 pc_info = &c.get_class_info();
161 // Method still not found. This shouldn't happen because basic (the
162 // base class of the algebraic hierarchy) registers a method for
163 // print_context (the base class of the print context hierarchy),
164 // so if we end up here, there's something wrong with the class
166 throw (std::runtime_error(std::string("basic::print(): method for ") + class_name() + "/" + c.class_name() + " not found"));
171 pdt[id](*this, c, level);
175 /** Default output to stream. */
176 void basic::do_print(const print_context & c, unsigned level) const
178 c.s << "[" << class_name() << " object]";
181 /** Tree output to stream. */
182 void basic::do_print_tree(const print_tree & c, unsigned level) const
184 c.s << std::string(level, ' ') << class_name() << " @" << this
185 << std::hex << ", hash=0x" << hashvalue << ", flags=0x" << flags << std::dec;
187 c.s << ", nops=" << nops();
189 for (size_t i=0; i<nops(); ++i)
190 op(i).print(c, level + c.delta_indent);
193 /** Python parsable output to stream. */
194 void basic::do_print_python_repr(const print_python_repr & c, unsigned level) const
196 c.s << class_name() << "()";
199 /** Little wrapper around print to be called within a debugger.
200 * This is needed because you cannot call foo.print(cout) from within the
201 * debugger because it might not know what cout is. This method can be
202 * invoked with no argument and it will simply print to stdout.
205 * @see basic::dbgprinttree */
206 void basic::dbgprint() const
208 this->print(print_dflt(std::cerr));
209 std::cerr << std::endl;
212 /** Little wrapper around printtree to be called within a debugger.
214 * @see basic::dbgprint */
215 void basic::dbgprinttree() const
217 this->print(print_tree(std::cerr));
220 /** Return relative operator precedence (for parenthezing output). */
221 unsigned basic::precedence() const
226 /** Information about the object.
228 * @see class info_flags */
229 bool basic::info(unsigned inf) const
231 // all possible properties are false for basic objects
235 /** Number of operands/members. */
236 size_t basic::nops() const
238 // iterating from 0 to nops() on atomic objects should be an empty loop,
239 // and accessing their elements is a range error. Container objects should
244 /** Return operand/member at position i. */
245 ex basic::op(size_t i) const
247 throw(std::range_error(std::string("basic::op(): ") + class_name() + std::string(" has no operands")));
250 /** Return modifyable operand/member at position i. */
251 ex & basic::let_op(size_t i)
253 ensure_if_modifiable();
254 throw(std::range_error(std::string("basic::let_op(): ") + class_name() + std::string(" has no operands")));
257 ex basic::operator[](const ex & index) const
259 if (is_exactly_a<numeric>(index))
260 return op(static_cast<size_t>(ex_to<numeric>(index).to_int()));
262 throw(std::invalid_argument(std::string("non-numeric indices not supported by ") + class_name()));
265 ex basic::operator[](size_t i) const
270 ex & basic::operator[](const ex & index)
272 if (is_exactly_a<numeric>(index))
273 return let_op(ex_to<numeric>(index).to_int());
275 throw(std::invalid_argument(std::string("non-numeric indices not supported by ") + class_name()));
278 ex & basic::operator[](size_t i)
283 /** Test for occurrence of a pattern. An object 'has' a pattern if it matches
284 * the pattern itself or one of the children 'has' it. As a consequence
285 * (according to the definition of children) given e=x+y+z, e.has(x) is true
286 * but e.has(x+y) is false. */
287 bool basic::has(const ex & pattern, unsigned options) const
290 if (match(pattern, repl_lst))
292 for (size_t i=0; i<nops(); i++)
293 if (op(i).has(pattern, options))
299 /** Construct new expression by applying the specified function to all
300 * sub-expressions (one level only, not recursively). */
301 ex basic::map(map_function & f) const
308 for (size_t i=0; i<num; i++) {
309 const ex & o = op(i);
311 if (!are_ex_trivially_equal(o, n)) {
319 copy->setflag(status_flags::dynallocated);
320 copy->clearflag(status_flags::hash_calculated | status_flags::expanded);
326 /** Check whether this is a polynomial in the given variables. */
327 bool basic::is_polynomial(const ex & var) const
329 return !has(var) || is_equal(ex_to<basic>(var));
332 /** Return degree of highest power in object s. */
333 int basic::degree(const ex & s) const
335 return is_equal(ex_to<basic>(s)) ? 1 : 0;
338 /** Return degree of lowest power in object s. */
339 int basic::ldegree(const ex & s) const
341 return is_equal(ex_to<basic>(s)) ? 1 : 0;
344 /** Return coefficient of degree n in object s. */
345 ex basic::coeff(const ex & s, int n) const
347 if (is_equal(ex_to<basic>(s)))
348 return n==1 ? _ex1 : _ex0;
350 return n==0 ? *this : _ex0;
353 /** Sort expanded expression in terms of powers of some object(s).
354 * @param s object(s) to sort in
355 * @param distributed recursive or distributed form (only used when s is a list) */
356 ex basic::collect(const ex & s, bool distributed) const
361 // List of objects specified
365 return collect(s.op(0));
367 else if (distributed) {
372 const lst& l(ex_to<lst>(s));
376 for (const_iterator xi=x.begin(); xi!=x.end(); ++xi) {
379 for (lst::const_iterator li=l.begin(); li!=l.end(); ++li) {
380 int cexp = pre_coeff.degree(*li);
381 pre_coeff = pre_coeff.coeff(*li, cexp);
382 key *= pow(*li, cexp);
384 exmap::iterator ci = cmap.find(key);
385 if (ci != cmap.end())
386 ci->second += pre_coeff;
388 cmap.insert(exmap::value_type(key, pre_coeff));
392 for (exmap::const_iterator mi=cmap.begin(); mi != cmap.end(); ++mi)
393 resv.push_back((mi->first)*(mi->second));
394 return (new add(resv))->setflag(status_flags::dynallocated);
400 size_t n = s.nops() - 1;
411 // Only one object specified
412 for (int n=this->ldegree(s); n<=this->degree(s); ++n)
413 x += this->coeff(s,n)*power(s,n);
416 // correct for lost fractional arguments and return
417 return x + (*this - x).expand();
420 /** Perform automatic non-interruptive term rewriting rules. */
421 ex basic::eval(int level) const
423 // There is nothing to do for basic objects:
427 /** Function object to be applied by basic::evalf(). */
428 struct evalf_map_function : public map_function {
430 evalf_map_function(int l) : level(l) {}
431 ex operator()(const ex & e) { return evalf(e, level); }
434 /** Evaluate object numerically. */
435 ex basic::evalf(int level) const
442 else if (level == -max_recursion_level)
443 throw(std::runtime_error("max recursion level reached"));
445 evalf_map_function map_evalf(level - 1);
446 return map(map_evalf);
451 /** Function object to be applied by basic::evalm(). */
452 struct evalm_map_function : public map_function {
453 ex operator()(const ex & e) { return evalm(e); }
456 /** Evaluate sums, products and integer powers of matrices. */
457 ex basic::evalm() const
462 return map(map_evalm);
465 /** Function object to be applied by basic::eval_integ(). */
466 struct eval_integ_map_function : public map_function {
467 ex operator()(const ex & e) { return eval_integ(e); }
470 /** Evaluate integrals, if result is known. */
471 ex basic::eval_integ() const
476 return map(map_eval_integ);
479 /** Perform automatic symbolic evaluations on indexed expression that
480 * contains this object as the base expression. */
481 ex basic::eval_indexed(const basic & i) const
482 // this function can't take a "const ex & i" because that would result
483 // in an infinite eval() loop
485 // There is nothing to do for basic objects
489 /** Add two indexed expressions. They are guaranteed to be of class indexed
490 * (or a subclass) and their indices are compatible. This function is used
491 * internally by simplify_indexed().
493 * @param self First indexed expression; its base object is *this
494 * @param other Second indexed expression
495 * @return sum of self and other
496 * @see ex::simplify_indexed() */
497 ex basic::add_indexed(const ex & self, const ex & other) const
502 /** Multiply an indexed expression with a scalar. This function is used
503 * internally by simplify_indexed().
505 * @param self Indexed expression; its base object is *this
506 * @param other Numeric value
507 * @return product of self and other
508 * @see ex::simplify_indexed() */
509 ex basic::scalar_mul_indexed(const ex & self, const numeric & other) const
514 /** Try to contract two indexed expressions that appear in the same product.
515 * If a contraction exists, the function overwrites one or both of the
516 * expressions and returns true. Otherwise it returns false. It is
517 * guaranteed that both expressions are of class indexed (or a subclass)
518 * and that at least one dummy index has been found. This functions is
519 * used internally by simplify_indexed().
521 * @param self Pointer to first indexed expression; its base object is *this
522 * @param other Pointer to second indexed expression
523 * @param v The complete vector of factors
524 * @return true if the contraction was successful, false otherwise
525 * @see ex::simplify_indexed() */
526 bool basic::contract_with(exvector::iterator self, exvector::iterator other, exvector & v) const
532 /** Check whether the expression matches a given pattern. For every wildcard
533 * object in the pattern, a pair with the wildcard as a key and matching
534 * expression as a value is added to repl_lst. */
535 bool basic::match(const ex & pattern, exmap& repl_lst) const
538 Sweet sweet shapes, sweet sweet shapes,
539 That's the key thing, right right.
540 Feed feed face, feed feed shapes,
541 But who is the king tonight?
542 Who is the king tonight?
543 Pattern is the thing, the key thing-a-ling,
544 But who is the king of Pattern?
545 But who is the king, the king thing-a-ling,
546 Who is the king of Pattern?
547 Bog is the king, the king thing-a-ling,
548 Bog is the king of Pattern.
549 Ba bu-bu-bu-bu bu-bu-bu-bu-bu-bu bu-bu
550 Bog is the king of Pattern.
553 if (is_exactly_a<wildcard>(pattern)) {
555 // Wildcard matches anything, but check whether we already have found
556 // a match for that wildcard first (if so, the earlier match must be
557 // the same expression)
558 for (exmap::const_iterator it = repl_lst.begin(); it != repl_lst.end(); ++it) {
559 if (it->first.is_equal(pattern))
560 return is_equal(ex_to<basic>(it->second));
562 repl_lst[pattern] = *this;
567 // Expression must be of the same type as the pattern
568 if (typeid(*this) != typeid(ex_to<basic>(pattern)))
571 // Number of subexpressions must match
572 if (nops() != pattern.nops())
575 // No subexpressions? Then just compare the objects (there can't be
576 // wildcards in the pattern)
578 return is_equal_same_type(ex_to<basic>(pattern));
580 // Check whether attributes that are not subexpressions match
581 if (!match_same_type(ex_to<basic>(pattern)))
584 // Even if the expression does not match the pattern, some of
585 // its subexpressions could match it. For example, x^5*y^(-1)
586 // does not match the pattern $0^5, but its subexpression x^5
587 // does. So, save repl_lst in order to not add bogus entries.
588 exmap tmp_repl = repl_lst;
589 // Otherwise the subexpressions must match one-to-one
590 for (size_t i=0; i<nops(); i++)
591 if (!op(i).match(pattern.op(i), tmp_repl))
594 // Looks similar enough, match found
600 /** Helper function for subs(). Does not recurse into subexpressions. */
601 ex basic::subs_one_level(const exmap & m, unsigned options) const
603 exmap::const_iterator it;
605 if (options & subs_options::no_pattern) {
612 for (it = m.begin(); it != m.end(); ++it) {
614 if (match(ex_to<basic>(it->first), repl_lst))
615 return it->second.subs(repl_lst, options | subs_options::no_pattern);
616 // avoid infinite recursion when re-substituting the wildcards
623 /** Substitute a set of objects by arbitrary expressions. The ex returned
624 * will already be evaluated. */
625 ex basic::subs(const exmap & m, unsigned options) const
630 // Substitute in subexpressions
631 for (size_t i=0; i<num; i++) {
632 const ex & orig_op = op(i);
633 const ex & subsed_op = orig_op.subs(m, options);
634 if (!are_ex_trivially_equal(orig_op, subsed_op)) {
636 // Something changed, clone the object
637 basic *copy = duplicate();
638 copy->setflag(status_flags::dynallocated);
639 copy->clearflag(status_flags::hash_calculated | status_flags::expanded);
641 // Substitute the changed operand
642 copy->let_op(i++) = subsed_op;
644 // Substitute the other operands
646 copy->let_op(i) = op(i).subs(m, options);
648 // Perform substitutions on the new object as a whole
649 return copy->subs_one_level(m, options);
654 // Nothing changed or no subexpressions
655 return subs_one_level(m, options);
658 /** Default interface of nth derivative ex::diff(s, n). It should be called
659 * instead of ::derivative(s) for first derivatives and for nth derivatives it
660 * just recurses down.
662 * @param s symbol to differentiate in
663 * @param nth order of differentiation
665 ex basic::diff(const symbol & s, unsigned nth) const
667 // trivial: zeroth derivative
671 // evaluate unevaluated *this before differentiating
672 if (!(flags & status_flags::evaluated))
673 return ex(*this).diff(s, nth);
675 ex ndiff = this->derivative(s);
676 while (!ndiff.is_zero() && // stop differentiating zeros
678 ndiff = ndiff.diff(s);
684 /** Return a vector containing the free indices of an expression. */
685 exvector basic::get_free_indices() const
687 return exvector(); // return an empty exvector
690 ex basic::conjugate() const
695 ex basic::real_part() const
697 return real_part_function(*this).hold();
700 ex basic::imag_part() const
702 return imag_part_function(*this).hold();
705 ex basic::eval_ncmul(const exvector & v) const
707 return hold_ncmul(v);
712 /** Function object to be applied by basic::derivative(). */
713 struct derivative_map_function : public map_function {
715 derivative_map_function(const symbol &sym) : s(sym) {}
716 ex operator()(const ex & e) { return diff(e, s); }
719 /** Default implementation of ex::diff(). It maps the operation on the
720 * operands (or returns 0 when the object has no operands).
723 ex basic::derivative(const symbol & s) const
728 derivative_map_function map_derivative(s);
729 return map(map_derivative);
733 /** Returns order relation between two objects of same type. This needs to be
734 * implemented by each class. It may never return anything else than 0,
735 * signalling equality, or +1 and -1 signalling inequality and determining
736 * the canonical ordering. (Perl hackers will wonder why C++ doesn't feature
737 * the spaceship operator <=> for denoting just this.) */
738 int basic::compare_same_type(const basic & other) const
740 return compare_pointers(this, &other);
743 /** Returns true if two objects of same type are equal. Normally needs
744 * not be reimplemented as long as it wasn't overwritten by some parent
745 * class, since it just calls compare_same_type(). The reason why this
746 * function exists is that sometimes it is easier to determine equality
747 * than an order relation and then it can be overridden. */
748 bool basic::is_equal_same_type(const basic & other) const
750 return compare_same_type(other)==0;
753 /** Returns true if the attributes of two objects are similar enough for
754 * a match. This function must not match subexpressions (this is already
755 * done by basic::match()). Only attributes not accessible by op() should
756 * be compared. This is also the reason why this function doesn't take the
757 * wildcard replacement list from match() as an argument: only subexpressions
758 * are subject to wildcard matches. Also, this function only needs to be
759 * implemented for container classes because is_equal_same_type() is
760 * automatically used instead of match_same_type() if nops() == 0.
762 * @see basic::match */
763 bool basic::match_same_type(const basic & other) const
765 // The default is to only consider subexpressions, but not any other
770 unsigned basic::return_type() const
772 return return_types::commutative;
775 return_type_t basic::return_type_tinfo() const
778 rt.tinfo = &typeid(*this);
783 /** Compute the hash value of an object and if it makes sense to store it in
784 * the objects status_flags, do so. The method inherited from class basic
785 * computes a hash value based on the type and hash values of possible
786 * members. For this reason it is well suited for container classes but
787 * atomic classes should override this implementation because otherwise they
788 * would all end up with the same hashvalue. */
789 unsigned basic::calchash() const
791 const void* this_tinfo = (const void*)typeid(*this).name();
792 unsigned v = golden_ratio_hash((p_int)this_tinfo);
793 for (size_t i=0; i<nops(); i++) {
795 v ^= this->op(i).gethash();
798 // store calculated hash value only if object is already evaluated
799 if (flags & status_flags::evaluated) {
800 setflag(status_flags::hash_calculated);
807 /** Function object to be applied by basic::expand(). */
808 struct expand_map_function : public map_function {
810 expand_map_function(unsigned o) : options(o) {}
811 ex operator()(const ex & e) { return e.expand(options); }
814 /** Expand expression, i.e. multiply it out and return the result as a new
816 ex basic::expand(unsigned options) const
819 return (options == 0) ? setflag(status_flags::expanded) : *this;
821 expand_map_function map_expand(options);
822 return ex_to<basic>(map(map_expand)).setflag(options == 0 ? status_flags::expanded : 0);
828 // non-virtual functions in this class
833 /** Compare objects syntactically to establish canonical ordering.
834 * All compare functions return: -1 for *this less than other, 0 equal,
836 int basic::compare(const basic & other) const
838 #ifdef GINAC_COMPARE_STATISTICS
839 compare_statistics.total_basic_compares++;
841 const unsigned hash_this = gethash();
842 const unsigned hash_other = other.gethash();
843 if (hash_this<hash_other) return -1;
844 if (hash_this>hash_other) return 1;
845 #ifdef GINAC_COMPARE_STATISTICS
846 compare_statistics.compare_same_hashvalue++;
849 const std::type_info& typeid_this = typeid(*this);
850 const std::type_info& typeid_other = typeid(other);
851 if (typeid_this == typeid_other) {
852 // int cmpval = compare_same_type(other);
854 // std::cout << "hash collision, same type: "
855 // << *this << " and " << other << std::endl;
856 // this->print(print_tree(std::cout));
857 // std::cout << " and ";
858 // other.print(print_tree(std::cout));
859 // std::cout << std::endl;
862 #ifdef GINAC_COMPARE_STATISTICS
863 compare_statistics.compare_same_type++;
865 return compare_same_type(other);
867 // std::cout << "hash collision, different types: "
868 // << *this << " and " << other << std::endl;
869 // this->print(print_tree(std::cout));
870 // std::cout << " and ";
871 // other.print(print_tree(std::cout));
872 // std::cout << std::endl;
873 return (typeid_this.before(typeid_other) ? -1 : 1);
877 /** Test for syntactic equality.
878 * This is only a quick test, meaning objects should be in the same domain.
879 * You might have to .expand(), .normal() objects first, depending on the
880 * domain of your computation, to get a more reliable answer.
882 * @see is_equal_same_type */
883 bool basic::is_equal(const basic & other) const
885 #ifdef GINAC_COMPARE_STATISTICS
886 compare_statistics.total_basic_is_equals++;
888 if (this->gethash()!=other.gethash())
890 #ifdef GINAC_COMPARE_STATISTICS
891 compare_statistics.is_equal_same_hashvalue++;
893 if (typeid(*this) != typeid(other))
896 #ifdef GINAC_COMPARE_STATISTICS
897 compare_statistics.is_equal_same_type++;
899 return is_equal_same_type(other);
904 /** Stop further evaluation.
906 * @see basic::eval */
907 const basic & basic::hold() const
909 return setflag(status_flags::evaluated);
912 /** Ensure the object may be modified without hurting others, throws if this
913 * is not the case. */
914 void basic::ensure_if_modifiable() const
916 if (get_refcount() > 1)
917 throw(std::runtime_error("cannot modify multiply referenced object"));
918 clearflag(status_flags::hash_calculated | status_flags::evaluated);
925 int max_recursion_level = 1024;
928 #ifdef GINAC_COMPARE_STATISTICS
929 compare_statistics_t::~compare_statistics_t()
931 std::clog << "ex::compare() called " << total_compares << " times" << std::endl;
932 std::clog << "nontrivial compares: " << nontrivial_compares << " times" << std::endl;
933 std::clog << "basic::compare() called " << total_basic_compares << " times" << std::endl;
934 std::clog << "same hashvalue in compare(): " << compare_same_hashvalue << " times" << std::endl;
935 std::clog << "compare_same_type() called " << compare_same_type << " times" << std::endl;
936 std::clog << std::endl;
937 std::clog << "ex::is_equal() called " << total_is_equals << " times" << std::endl;
938 std::clog << "nontrivial is_equals: " << nontrivial_is_equals << " times" << std::endl;
939 std::clog << "basic::is_equal() called " << total_basic_is_equals << " times" << std::endl;
940 std::clog << "same hashvalue in is_equal(): " << is_equal_same_hashvalue << " times" << std::endl;
941 std::clog << "is_equal_same_type() called " << is_equal_same_type << " times" << std::endl;
942 std::clog << std::endl;
943 std::clog << "basic::gethash() called " << total_gethash << " times" << std::endl;
944 std::clog << "used cached hashvalue " << gethash_cached << " times" << std::endl;
947 compare_statistics_t compare_statistics;