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 void basic::read_archive(const archive_node& n, lst& syms)
98 /** Archive the object. */
99 void basic::archive(archive_node &n) const
101 n.add_string("class", class_name());
105 // new virtual functions which can be overridden by derived classes
110 /** Output to stream. This performs double dispatch on the dynamic type of
111 * *this and the dynamic type of the supplied print context.
112 * @param c print context object that describes the output formatting
113 * @param level value that is used to identify the precedence or indentation
114 * level for placing parentheses and formatting */
115 void basic::print(const print_context & c, unsigned level) const
117 print_dispatch(get_class_info(), c, level);
120 /** Like print(), but dispatch to the specified class. Can be used by
121 * implementations of print methods to dispatch to the method of the
124 * @see basic::print */
125 void basic::print_dispatch(const registered_class_info & ri, const print_context & c, unsigned level) const
127 // Double dispatch on object type and print_context type
128 const registered_class_info * reg_info = &ri;
129 const print_context_class_info * pc_info = &c.get_class_info();
132 const std::vector<print_functor> & pdt = reg_info->options.get_print_dispatch_table();
135 unsigned id = pc_info->options.get_id();
136 if (id >= pdt.size() || !(pdt[id].is_valid())) {
138 // Method not found, try parent print_context class
139 const print_context_class_info * parent_pc_info = pc_info->get_parent();
140 if (parent_pc_info) {
141 pc_info = parent_pc_info;
145 // Method still not found, try parent class
146 const registered_class_info * parent_reg_info = reg_info->get_parent();
147 if (parent_reg_info) {
148 reg_info = parent_reg_info;
149 pc_info = &c.get_class_info();
153 // Method still not found. This shouldn't happen because basic (the
154 // base class of the algebraic hierarchy) registers a method for
155 // print_context (the base class of the print context hierarchy),
156 // so if we end up here, there's something wrong with the class
158 throw (std::runtime_error(std::string("basic::print(): method for ") + class_name() + "/" + c.class_name() + " not found"));
163 pdt[id](*this, c, level);
167 /** Default output to stream. */
168 void basic::do_print(const print_context & c, unsigned level) const
170 c.s << "[" << class_name() << " object]";
173 /** Tree output to stream. */
174 void basic::do_print_tree(const print_tree & c, unsigned level) const
176 c.s << std::string(level, ' ') << class_name() << " @" << this
177 << std::hex << ", hash=0x" << hashvalue << ", flags=0x" << flags << std::dec;
179 c.s << ", nops=" << nops();
181 for (size_t i=0; i<nops(); ++i)
182 op(i).print(c, level + c.delta_indent);
185 /** Python parsable output to stream. */
186 void basic::do_print_python_repr(const print_python_repr & c, unsigned level) const
188 c.s << class_name() << "()";
191 /** Little wrapper around print to be called within a debugger.
192 * This is needed because you cannot call foo.print(cout) from within the
193 * debugger because it might not know what cout is. This method can be
194 * invoked with no argument and it will simply print to stdout.
197 * @see basic::dbgprinttree */
198 void basic::dbgprint() const
200 this->print(print_dflt(std::cerr));
201 std::cerr << std::endl;
204 /** Little wrapper around printtree to be called within a debugger.
206 * @see basic::dbgprint */
207 void basic::dbgprinttree() const
209 this->print(print_tree(std::cerr));
212 /** Return relative operator precedence (for parenthezing output). */
213 unsigned basic::precedence() const
218 /** Information about the object.
220 * @see class info_flags */
221 bool basic::info(unsigned inf) const
223 // all possible properties are false for basic objects
227 /** Number of operands/members. */
228 size_t basic::nops() const
230 // iterating from 0 to nops() on atomic objects should be an empty loop,
231 // and accessing their elements is a range error. Container objects should
236 /** Return operand/member at position i. */
237 ex basic::op(size_t i) const
239 throw(std::range_error(std::string("basic::op(): ") + class_name() + std::string(" has no operands")));
242 /** Return modifyable operand/member at position i. */
243 ex & basic::let_op(size_t i)
245 ensure_if_modifiable();
246 throw(std::range_error(std::string("basic::let_op(): ") + class_name() + std::string(" has no operands")));
249 ex basic::operator[](const ex & index) const
251 if (is_exactly_a<numeric>(index))
252 return op(static_cast<size_t>(ex_to<numeric>(index).to_int()));
254 throw(std::invalid_argument(std::string("non-numeric indices not supported by ") + class_name()));
257 ex basic::operator[](size_t i) const
262 ex & basic::operator[](const ex & index)
264 if (is_exactly_a<numeric>(index))
265 return let_op(ex_to<numeric>(index).to_int());
267 throw(std::invalid_argument(std::string("non-numeric indices not supported by ") + class_name()));
270 ex & basic::operator[](size_t i)
275 /** Test for occurrence of a pattern. An object 'has' a pattern if it matches
276 * the pattern itself or one of the children 'has' it. As a consequence
277 * (according to the definition of children) given e=x+y+z, e.has(x) is true
278 * but e.has(x+y) is false. */
279 bool basic::has(const ex & pattern, unsigned options) const
282 if (match(pattern, repl_lst))
284 for (size_t i=0; i<nops(); i++)
285 if (op(i).has(pattern, options))
291 /** Construct new expression by applying the specified function to all
292 * sub-expressions (one level only, not recursively). */
293 ex basic::map(map_function & f) const
300 for (size_t i=0; i<num; i++) {
301 const ex & o = op(i);
303 if (!are_ex_trivially_equal(o, n)) {
311 copy->setflag(status_flags::dynallocated);
312 copy->clearflag(status_flags::hash_calculated | status_flags::expanded);
318 /** Check whether this is a polynomial in the given variables. */
319 bool basic::is_polynomial(const ex & var) const
321 return !has(var) || is_equal(ex_to<basic>(var));
324 /** Return degree of highest power in object s. */
325 int basic::degree(const ex & s) const
327 return is_equal(ex_to<basic>(s)) ? 1 : 0;
330 /** Return degree of lowest power in object s. */
331 int basic::ldegree(const ex & s) const
333 return is_equal(ex_to<basic>(s)) ? 1 : 0;
336 /** Return coefficient of degree n in object s. */
337 ex basic::coeff(const ex & s, int n) const
339 if (is_equal(ex_to<basic>(s)))
340 return n==1 ? _ex1 : _ex0;
342 return n==0 ? *this : _ex0;
345 /** Sort expanded expression in terms of powers of some object(s).
346 * @param s object(s) to sort in
347 * @param distributed recursive or distributed form (only used when s is a list) */
348 ex basic::collect(const ex & s, bool distributed) const
353 // List of objects specified
357 return collect(s.op(0));
359 else if (distributed) {
364 const lst& l(ex_to<lst>(s));
368 for (const_iterator xi=x.begin(); xi!=x.end(); ++xi) {
371 for (lst::const_iterator li=l.begin(); li!=l.end(); ++li) {
372 int cexp = pre_coeff.degree(*li);
373 pre_coeff = pre_coeff.coeff(*li, cexp);
374 key *= pow(*li, cexp);
376 exmap::iterator ci = cmap.find(key);
377 if (ci != cmap.end())
378 ci->second += pre_coeff;
380 cmap.insert(exmap::value_type(key, pre_coeff));
384 for (exmap::const_iterator mi=cmap.begin(); mi != cmap.end(); ++mi)
385 resv.push_back((mi->first)*(mi->second));
386 return (new add(resv))->setflag(status_flags::dynallocated);
392 size_t n = s.nops() - 1;
403 // Only one object specified
404 for (int n=this->ldegree(s); n<=this->degree(s); ++n)
405 x += this->coeff(s,n)*power(s,n);
408 // correct for lost fractional arguments and return
409 return x + (*this - x).expand();
412 /** Perform automatic non-interruptive term rewriting rules. */
413 ex basic::eval(int level) const
415 // There is nothing to do for basic objects:
419 /** Function object to be applied by basic::evalf(). */
420 struct evalf_map_function : public map_function {
422 evalf_map_function(int l) : level(l) {}
423 ex operator()(const ex & e) { return evalf(e, level); }
426 /** Evaluate object numerically. */
427 ex basic::evalf(int level) const
434 else if (level == -max_recursion_level)
435 throw(std::runtime_error("max recursion level reached"));
437 evalf_map_function map_evalf(level - 1);
438 return map(map_evalf);
443 /** Function object to be applied by basic::evalm(). */
444 struct evalm_map_function : public map_function {
445 ex operator()(const ex & e) { return evalm(e); }
448 /** Evaluate sums, products and integer powers of matrices. */
449 ex basic::evalm() const
454 return map(map_evalm);
457 /** Function object to be applied by basic::eval_integ(). */
458 struct eval_integ_map_function : public map_function {
459 ex operator()(const ex & e) { return eval_integ(e); }
462 /** Evaluate integrals, if result is known. */
463 ex basic::eval_integ() const
468 return map(map_eval_integ);
471 /** Perform automatic symbolic evaluations on indexed expression that
472 * contains this object as the base expression. */
473 ex basic::eval_indexed(const basic & i) const
474 // this function can't take a "const ex & i" because that would result
475 // in an infinite eval() loop
477 // There is nothing to do for basic objects
481 /** Add two indexed expressions. They are guaranteed to be of class indexed
482 * (or a subclass) and their indices are compatible. This function is used
483 * internally by simplify_indexed().
485 * @param self First indexed expression; its base object is *this
486 * @param other Second indexed expression
487 * @return sum of self and other
488 * @see ex::simplify_indexed() */
489 ex basic::add_indexed(const ex & self, const ex & other) const
494 /** Multiply an indexed expression with a scalar. This function is used
495 * internally by simplify_indexed().
497 * @param self Indexed expression; its base object is *this
498 * @param other Numeric value
499 * @return product of self and other
500 * @see ex::simplify_indexed() */
501 ex basic::scalar_mul_indexed(const ex & self, const numeric & other) const
506 /** Try to contract two indexed expressions that appear in the same product.
507 * If a contraction exists, the function overwrites one or both of the
508 * expressions and returns true. Otherwise it returns false. It is
509 * guaranteed that both expressions are of class indexed (or a subclass)
510 * and that at least one dummy index has been found. This functions is
511 * used internally by simplify_indexed().
513 * @param self Pointer to first indexed expression; its base object is *this
514 * @param other Pointer to second indexed expression
515 * @param v The complete vector of factors
516 * @return true if the contraction was successful, false otherwise
517 * @see ex::simplify_indexed() */
518 bool basic::contract_with(exvector::iterator self, exvector::iterator other, exvector & v) const
524 /** Check whether the expression matches a given pattern. For every wildcard
525 * object in the pattern, a pair with the wildcard as a key and matching
526 * expression as a value is added to repl_lst. */
527 bool basic::match(const ex & pattern, exmap& repl_lst) const
530 Sweet sweet shapes, sweet sweet shapes,
531 That's the key thing, right right.
532 Feed feed face, feed feed shapes,
533 But who is the king tonight?
534 Who is the king tonight?
535 Pattern is the thing, the key thing-a-ling,
536 But who is the king of Pattern?
537 But who is the king, the king thing-a-ling,
538 Who is the king of Pattern?
539 Bog is the king, the king thing-a-ling,
540 Bog is the king of Pattern.
541 Ba bu-bu-bu-bu bu-bu-bu-bu-bu-bu bu-bu
542 Bog is the king of Pattern.
545 if (is_exactly_a<wildcard>(pattern)) {
547 // Wildcard matches anything, but check whether we already have found
548 // a match for that wildcard first (if so, the earlier match must be
549 // the same expression)
550 for (exmap::const_iterator it = repl_lst.begin(); it != repl_lst.end(); ++it) {
551 if (it->first.is_equal(pattern))
552 return is_equal(ex_to<basic>(it->second));
554 repl_lst[pattern] = *this;
559 // Expression must be of the same type as the pattern
560 if (typeid(*this) != typeid(ex_to<basic>(pattern)))
563 // Number of subexpressions must match
564 if (nops() != pattern.nops())
567 // No subexpressions? Then just compare the objects (there can't be
568 // wildcards in the pattern)
570 return is_equal_same_type(ex_to<basic>(pattern));
572 // Check whether attributes that are not subexpressions match
573 if (!match_same_type(ex_to<basic>(pattern)))
576 // Even if the expression does not match the pattern, some of
577 // its subexpressions could match it. For example, x^5*y^(-1)
578 // does not match the pattern $0^5, but its subexpression x^5
579 // does. So, save repl_lst in order to not add bogus entries.
580 exmap tmp_repl = repl_lst;
581 // Otherwise the subexpressions must match one-to-one
582 for (size_t i=0; i<nops(); i++)
583 if (!op(i).match(pattern.op(i), tmp_repl))
586 // Looks similar enough, match found
592 /** Helper function for subs(). Does not recurse into subexpressions. */
593 ex basic::subs_one_level(const exmap & m, unsigned options) const
595 exmap::const_iterator it;
597 if (options & subs_options::no_pattern) {
604 for (it = m.begin(); it != m.end(); ++it) {
606 if (match(ex_to<basic>(it->first), repl_lst))
607 return it->second.subs(repl_lst, options | subs_options::no_pattern);
608 // avoid infinite recursion when re-substituting the wildcards
615 /** Substitute a set of objects by arbitrary expressions. The ex returned
616 * will already be evaluated. */
617 ex basic::subs(const exmap & m, unsigned options) const
622 // Substitute in subexpressions
623 for (size_t i=0; i<num; i++) {
624 const ex & orig_op = op(i);
625 const ex & subsed_op = orig_op.subs(m, options);
626 if (!are_ex_trivially_equal(orig_op, subsed_op)) {
628 // Something changed, clone the object
629 basic *copy = duplicate();
630 copy->setflag(status_flags::dynallocated);
631 copy->clearflag(status_flags::hash_calculated | status_flags::expanded);
633 // Substitute the changed operand
634 copy->let_op(i++) = subsed_op;
636 // Substitute the other operands
638 copy->let_op(i) = op(i).subs(m, options);
640 // Perform substitutions on the new object as a whole
641 return copy->subs_one_level(m, options);
646 // Nothing changed or no subexpressions
647 return subs_one_level(m, options);
650 /** Default interface of nth derivative ex::diff(s, n). It should be called
651 * instead of ::derivative(s) for first derivatives and for nth derivatives it
652 * just recurses down.
654 * @param s symbol to differentiate in
655 * @param nth order of differentiation
657 ex basic::diff(const symbol & s, unsigned nth) const
659 // trivial: zeroth derivative
663 // evaluate unevaluated *this before differentiating
664 if (!(flags & status_flags::evaluated))
665 return ex(*this).diff(s, nth);
667 ex ndiff = this->derivative(s);
668 while (!ndiff.is_zero() && // stop differentiating zeros
670 ndiff = ndiff.diff(s);
676 /** Return a vector containing the free indices of an expression. */
677 exvector basic::get_free_indices() const
679 return exvector(); // return an empty exvector
682 ex basic::conjugate() const
687 ex basic::real_part() const
689 return real_part_function(*this).hold();
692 ex basic::imag_part() const
694 return imag_part_function(*this).hold();
697 ex basic::eval_ncmul(const exvector & v) const
699 return hold_ncmul(v);
704 /** Function object to be applied by basic::derivative(). */
705 struct derivative_map_function : public map_function {
707 derivative_map_function(const symbol &sym) : s(sym) {}
708 ex operator()(const ex & e) { return diff(e, s); }
711 /** Default implementation of ex::diff(). It maps the operation on the
712 * operands (or returns 0 when the object has no operands).
715 ex basic::derivative(const symbol & s) const
720 derivative_map_function map_derivative(s);
721 return map(map_derivative);
725 /** Returns order relation between two objects of same type. This needs to be
726 * implemented by each class. It may never return anything else than 0,
727 * signalling equality, or +1 and -1 signalling inequality and determining
728 * the canonical ordering. (Perl hackers will wonder why C++ doesn't feature
729 * the spaceship operator <=> for denoting just this.) */
730 int basic::compare_same_type(const basic & other) const
732 return compare_pointers(this, &other);
735 /** Returns true if two objects of same type are equal. Normally needs
736 * not be reimplemented as long as it wasn't overwritten by some parent
737 * class, since it just calls compare_same_type(). The reason why this
738 * function exists is that sometimes it is easier to determine equality
739 * than an order relation and then it can be overridden. */
740 bool basic::is_equal_same_type(const basic & other) const
742 return compare_same_type(other)==0;
745 /** Returns true if the attributes of two objects are similar enough for
746 * a match. This function must not match subexpressions (this is already
747 * done by basic::match()). Only attributes not accessible by op() should
748 * be compared. This is also the reason why this function doesn't take the
749 * wildcard replacement list from match() as an argument: only subexpressions
750 * are subject to wildcard matches. Also, this function only needs to be
751 * implemented for container classes because is_equal_same_type() is
752 * automatically used instead of match_same_type() if nops() == 0.
754 * @see basic::match */
755 bool basic::match_same_type(const basic & other) const
757 // The default is to only consider subexpressions, but not any other
762 unsigned basic::return_type() const
764 return return_types::commutative;
767 return_type_t basic::return_type_tinfo() const
770 rt.tinfo = &typeid(*this);
775 /** Compute the hash value of an object and if it makes sense to store it in
776 * the objects status_flags, do so. The method inherited from class basic
777 * computes a hash value based on the type and hash values of possible
778 * members. For this reason it is well suited for container classes but
779 * atomic classes should override this implementation because otherwise they
780 * would all end up with the same hashvalue. */
781 unsigned basic::calchash() const
783 const void* this_tinfo = (const void*)typeid(*this).name();
784 unsigned v = golden_ratio_hash((p_int)this_tinfo);
785 for (size_t i=0; i<nops(); i++) {
787 v ^= this->op(i).gethash();
790 // store calculated hash value only if object is already evaluated
791 if (flags & status_flags::evaluated) {
792 setflag(status_flags::hash_calculated);
799 /** Function object to be applied by basic::expand(). */
800 struct expand_map_function : public map_function {
802 expand_map_function(unsigned o) : options(o) {}
803 ex operator()(const ex & e) { return e.expand(options); }
806 /** Expand expression, i.e. multiply it out and return the result as a new
808 ex basic::expand(unsigned options) const
811 return (options == 0) ? setflag(status_flags::expanded) : *this;
813 expand_map_function map_expand(options);
814 return ex_to<basic>(map(map_expand)).setflag(options == 0 ? status_flags::expanded : 0);
820 // non-virtual functions in this class
825 /** Compare objects syntactically to establish canonical ordering.
826 * All compare functions return: -1 for *this less than other, 0 equal,
828 int basic::compare(const basic & other) const
830 #ifdef GINAC_COMPARE_STATISTICS
831 compare_statistics.total_basic_compares++;
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++;
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);
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;
854 #ifdef GINAC_COMPARE_STATISTICS
855 compare_statistics.compare_same_type++;
857 return compare_same_type(other);
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);
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.
874 * @see is_equal_same_type */
875 bool basic::is_equal(const basic & other) const
877 #ifdef GINAC_COMPARE_STATISTICS
878 compare_statistics.total_basic_is_equals++;
880 if (this->gethash()!=other.gethash())
882 #ifdef GINAC_COMPARE_STATISTICS
883 compare_statistics.is_equal_same_hashvalue++;
885 if (typeid(*this) != typeid(other))
888 #ifdef GINAC_COMPARE_STATISTICS
889 compare_statistics.is_equal_same_type++;
891 return is_equal_same_type(other);
896 /** Stop further evaluation.
898 * @see basic::eval */
899 const basic & basic::hold() const
901 return setflag(status_flags::evaluated);
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
908 if (get_refcount() > 1)
909 throw(std::runtime_error("cannot modify multiply referenced object"));
910 clearflag(status_flags::hash_calculated | status_flags::evaluated);
917 int max_recursion_level = 1024;
920 #ifdef GINAC_COMPARE_STATISTICS
921 compare_statistics_t::~compare_statistics_t()
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;
939 compare_statistics_t compare_statistics;