3 * Implementation of GiNaC's ABC. */
6 * GiNaC Copyright (C) 1999-2015 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
31 #include "relational.h"
32 #include "operators.h"
36 #include "hash_seed.h"
45 GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(basic, void,
46 print_func<print_context>(&basic::do_print).
47 print_func<print_tree>(&basic::do_print_tree).
48 print_func<print_python_repr>(&basic::do_print_python_repr))
51 // default constructor, destructor, copy constructor and assignment operator
56 /** basic copy constructor: implicitly assumes that the other class is of
57 * the exact same type (as it's used by duplicate()), so it can copy the
58 * tinfo_key and the hash value. */
59 basic::basic(const basic & other) : flags(other.flags & ~status_flags::dynallocated), hashvalue(other.hashvalue)
63 /** basic assignment operator: the other object might be of a derived class. */
64 const basic & basic::operator=(const basic & other)
66 unsigned fl = other.flags & ~status_flags::dynallocated;
67 if (typeid(*this) != typeid(other)) {
68 // The other object is of a derived class, so clear the flags as they
69 // might no longer apply (especially hash_calculated). Oh, and don't
70 // copy the tinfo_key: it is already set correctly for this object.
71 fl &= ~(status_flags::evaluated | status_flags::expanded | status_flags::hash_calculated);
73 // The objects are of the exact same class, so copy the hash value.
74 hashvalue = other.hashvalue;
95 /** Construct object from archive_node. */
96 void basic::read_archive(const archive_node& n, lst& syms)
99 /** Archive the object. */
100 void basic::archive(archive_node &n) const
102 n.add_string("class", class_name());
106 // new virtual functions which can be overridden by derived classes
111 /** Output to stream. This performs double dispatch on the dynamic type of
112 * *this and the dynamic type of the supplied print context.
113 * @param c print context object that describes the output formatting
114 * @param level value that is used to identify the precedence or indentation
115 * level for placing parentheses and formatting */
116 void basic::print(const print_context & c, unsigned level) const
118 print_dispatch(get_class_info(), c, level);
121 /** Like print(), but dispatch to the specified class. Can be used by
122 * implementations of print methods to dispatch to the method of the
125 * @see basic::print */
126 void basic::print_dispatch(const registered_class_info & ri, const print_context & c, unsigned level) const
128 // Double dispatch on object type and print_context type
129 const registered_class_info * reg_info = &ri;
130 const print_context_class_info * pc_info = &c.get_class_info();
133 const std::vector<print_functor> & pdt = reg_info->options.get_print_dispatch_table();
136 unsigned id = pc_info->options.get_id();
137 if (id >= pdt.size() || !(pdt[id].is_valid())) {
139 // Method not found, try parent print_context class
140 const print_context_class_info * parent_pc_info = pc_info->get_parent();
141 if (parent_pc_info) {
142 pc_info = parent_pc_info;
146 // Method still not found, try parent class
147 const registered_class_info * parent_reg_info = reg_info->get_parent();
148 if (parent_reg_info) {
149 reg_info = parent_reg_info;
150 pc_info = &c.get_class_info();
154 // Method still not found. This shouldn't happen because basic (the
155 // base class of the algebraic hierarchy) registers a method for
156 // print_context (the base class of the print context hierarchy),
157 // so if we end up here, there's something wrong with the class
159 throw (std::runtime_error(std::string("basic::print(): method for ") + class_name() + "/" + c.class_name() + " not found"));
164 pdt[id](*this, c, level);
168 /** Default output to stream. */
169 void basic::do_print(const print_context & c, unsigned level) const
171 c.s << "[" << class_name() << " object]";
174 /** Tree output to stream. */
175 void basic::do_print_tree(const print_tree & c, unsigned level) const
177 c.s << std::string(level, ' ') << class_name() << " @" << this
178 << std::hex << ", hash=0x" << hashvalue << ", flags=0x" << flags << std::dec;
180 c.s << ", nops=" << nops();
182 for (size_t i=0; i<nops(); ++i)
183 op(i).print(c, level + c.delta_indent);
186 /** Python parsable output to stream. */
187 void basic::do_print_python_repr(const print_python_repr & c, unsigned level) const
189 c.s << class_name() << "()";
192 /** Little wrapper around print to be called within a debugger.
193 * This is needed because you cannot call foo.print(cout) from within the
194 * debugger because it might not know what cout is. This method can be
195 * invoked with no argument and it will simply print to stdout.
198 * @see basic::dbgprinttree */
199 void basic::dbgprint() const
201 this->print(print_dflt(std::cerr));
202 std::cerr << std::endl;
205 /** Little wrapper around printtree to be called within a debugger.
207 * @see basic::dbgprint */
208 void basic::dbgprinttree() const
210 this->print(print_tree(std::cerr));
213 /** Return relative operator precedence (for parenthezing output). */
214 unsigned basic::precedence() const
219 /** Information about the object.
221 * @see class info_flags */
222 bool basic::info(unsigned inf) const
224 // all possible properties are false for basic objects
228 /** Number of operands/members. */
229 size_t basic::nops() const
231 // iterating from 0 to nops() on atomic objects should be an empty loop,
232 // and accessing their elements is a range error. Container objects should
237 /** Return operand/member at position i. */
238 ex basic::op(size_t i) const
240 throw(std::range_error(std::string("basic::op(): ") + class_name() + std::string(" has no operands")));
243 /** Return modifiable operand/member at position i. */
244 ex & basic::let_op(size_t i)
246 ensure_if_modifiable();
247 throw(std::range_error(std::string("basic::let_op(): ") + class_name() + std::string(" has no operands")));
250 ex basic::operator[](const ex & index) const
252 if (is_exactly_a<numeric>(index))
253 return op(static_cast<size_t>(ex_to<numeric>(index).to_int()));
255 throw(std::invalid_argument(std::string("non-numeric indices not supported by ") + class_name()));
258 ex basic::operator[](size_t i) const
263 ex & basic::operator[](const ex & index)
265 if (is_exactly_a<numeric>(index))
266 return let_op(ex_to<numeric>(index).to_int());
268 throw(std::invalid_argument(std::string("non-numeric indices not supported by ") + class_name()));
271 ex & basic::operator[](size_t i)
276 /** Test for occurrence of a pattern. An object 'has' a pattern if it matches
277 * the pattern itself or one of the children 'has' it. As a consequence
278 * (according to the definition of children) given e=x+y+z, e.has(x) is true
279 * but e.has(x+y) is false. */
280 bool basic::has(const ex & pattern, unsigned options) const
283 if (match(pattern, repl_lst))
285 for (size_t i=0; i<nops(); i++)
286 if (op(i).has(pattern, options))
292 /** Construct new expression by applying the specified function to all
293 * sub-expressions (one level only, not recursively). */
294 ex basic::map(map_function & f) const
300 basic *copy = nullptr;
301 for (size_t i=0; i<num; i++) {
302 const ex & o = op(i);
304 if (!are_ex_trivially_equal(o, n)) {
312 copy->setflag(status_flags::dynallocated);
313 copy->clearflag(status_flags::hash_calculated | status_flags::expanded);
319 /** Check whether this is a polynomial in the given variables. */
320 bool basic::is_polynomial(const ex & var) const
322 return !has(var) || is_equal(ex_to<basic>(var));
325 /** Return degree of highest power in object s. */
326 int basic::degree(const ex & s) const
328 return is_equal(ex_to<basic>(s)) ? 1 : 0;
331 /** Return degree of lowest power in object s. */
332 int basic::ldegree(const ex & s) const
334 return is_equal(ex_to<basic>(s)) ? 1 : 0;
337 /** Return coefficient of degree n in object s. */
338 ex basic::coeff(const ex & s, int n) const
340 if (is_equal(ex_to<basic>(s)))
341 return n==1 ? _ex1 : _ex0;
343 return n==0 ? *this : _ex0;
346 /** Sort expanded expression in terms of powers of some object(s).
347 * @param s object(s) to sort in
348 * @param distributed recursive or distributed form (only used when s is a list) */
349 ex basic::collect(const ex & s, bool distributed) const
354 // List of objects specified
358 return collect(s.op(0));
360 else if (distributed) {
365 const lst& l(ex_to<lst>(s));
369 for (const auto & xi : x) {
372 for (auto & li : l) {
373 int cexp = pre_coeff.degree(li);
374 pre_coeff = pre_coeff.coeff(li, cexp);
375 key *= pow(li, cexp);
377 exmap::iterator ci = cmap.find(key);
378 if (ci != cmap.end())
379 ci->second += pre_coeff;
381 cmap.insert(exmap::value_type(key, pre_coeff));
385 for (auto & mi : cmap)
386 resv.push_back((mi.first)*(mi.second));
387 return (new add(resv))->setflag(status_flags::dynallocated);
393 size_t n = s.nops() - 1;
404 // Only one object specified
405 for (int n=this->ldegree(s); n<=this->degree(s); ++n)
406 x += this->coeff(s,n)*power(s,n);
409 // correct for lost fractional arguments and return
410 return x + (*this - x).expand();
413 /** Perform automatic non-interruptive term rewriting rules. */
414 ex basic::eval(int level) const
416 // There is nothing to do for basic objects:
420 /** Function object to be applied by basic::evalf(). */
421 struct evalf_map_function : public map_function {
423 evalf_map_function(int l) : level(l) {}
424 ex operator()(const ex & e) { return evalf(e, level); }
427 /** Evaluate object numerically. */
428 ex basic::evalf(int level) const
435 else if (level == -max_recursion_level)
436 throw(std::runtime_error("max recursion level reached"));
438 evalf_map_function map_evalf(level - 1);
439 return map(map_evalf);
444 /** Function object to be applied by basic::evalm(). */
445 struct evalm_map_function : public map_function {
446 ex operator()(const ex & e) { return evalm(e); }
449 /** Evaluate sums, products and integer powers of matrices. */
450 ex basic::evalm() const
455 return map(map_evalm);
458 /** Function object to be applied by basic::eval_integ(). */
459 struct eval_integ_map_function : public map_function {
460 ex operator()(const ex & e) { return eval_integ(e); }
463 /** Evaluate integrals, if result is known. */
464 ex basic::eval_integ() const
469 return map(map_eval_integ);
472 /** Perform automatic symbolic evaluations on indexed expression that
473 * contains this object as the base expression. */
474 ex basic::eval_indexed(const basic & i) const
475 // this function can't take a "const ex & i" because that would result
476 // in an infinite eval() loop
478 // There is nothing to do for basic objects
482 /** Add two indexed expressions. They are guaranteed to be of class indexed
483 * (or a subclass) and their indices are compatible. This function is used
484 * internally by simplify_indexed().
486 * @param self First indexed expression; its base object is *this
487 * @param other Second indexed expression
488 * @return sum of self and other
489 * @see ex::simplify_indexed() */
490 ex basic::add_indexed(const ex & self, const ex & other) const
495 /** Multiply an indexed expression with a scalar. This function is used
496 * internally by simplify_indexed().
498 * @param self Indexed expression; its base object is *this
499 * @param other Numeric value
500 * @return product of self and other
501 * @see ex::simplify_indexed() */
502 ex basic::scalar_mul_indexed(const ex & self, const numeric & other) const
507 /** Try to contract two indexed expressions that appear in the same product.
508 * If a contraction exists, the function overwrites one or both of the
509 * expressions and returns true. Otherwise it returns false. It is
510 * guaranteed that both expressions are of class indexed (or a subclass)
511 * and that at least one dummy index has been found. This functions is
512 * used internally by simplify_indexed().
514 * @param self Pointer to first indexed expression; its base object is *this
515 * @param other Pointer to second indexed expression
516 * @param v The complete vector of factors
517 * @return true if the contraction was successful, false otherwise
518 * @see ex::simplify_indexed() */
519 bool basic::contract_with(exvector::iterator self, exvector::iterator other, exvector & v) const
525 /** Check whether the expression matches a given pattern. For every wildcard
526 * object in the pattern, a pair with the wildcard as a key and matching
527 * expression as a value is added to repl_lst. */
528 bool basic::match(const ex & pattern, exmap& repl_lst) const
531 Sweet sweet shapes, sweet sweet shapes,
532 That's the key thing, right right.
533 Feed feed face, feed feed shapes,
534 But who is the king tonight?
535 Who is the king tonight?
536 Pattern is the thing, the key thing-a-ling,
537 But who is the king of Pattern?
538 But who is the king, the king thing-a-ling,
539 Who is the king of Pattern?
540 Bog is the king, the king thing-a-ling,
541 Bog is the king of Pattern.
542 Ba bu-bu-bu-bu bu-bu-bu-bu-bu-bu bu-bu
543 Bog is the king of Pattern.
546 if (is_exactly_a<wildcard>(pattern)) {
548 // Wildcard matches anything, but check whether we already have found
549 // a match for that wildcard first (if so, the earlier match must be
550 // the same expression)
551 for (auto & it : repl_lst) {
552 if (it.first.is_equal(pattern))
553 return is_equal(ex_to<basic>(it.second));
555 repl_lst[pattern] = *this;
560 // Expression must be of the same type as the pattern
561 if (typeid(*this) != typeid(ex_to<basic>(pattern)))
564 // Number of subexpressions must match
565 if (nops() != pattern.nops())
568 // No subexpressions? Then just compare the objects (there can't be
569 // wildcards in the pattern)
571 return is_equal_same_type(ex_to<basic>(pattern));
573 // Check whether attributes that are not subexpressions match
574 if (!match_same_type(ex_to<basic>(pattern)))
577 // Even if the expression does not match the pattern, some of
578 // its subexpressions could match it. For example, x^5*y^(-1)
579 // does not match the pattern $0^5, but its subexpression x^5
580 // does. So, save repl_lst in order to not add bogus entries.
581 exmap tmp_repl = repl_lst;
582 // Otherwise the subexpressions must match one-to-one
583 for (size_t i=0; i<nops(); i++)
584 if (!op(i).match(pattern.op(i), tmp_repl))
587 // Looks similar enough, match found
593 /** Helper function for subs(). Does not recurse into subexpressions. */
594 ex basic::subs_one_level(const exmap & m, unsigned options) const
596 if (options & subs_options::no_pattern) {
597 auto it = m.find(*this);
602 for (auto & it : m) {
604 if (match(ex_to<basic>(it.first), repl_lst))
605 return it.second.subs(repl_lst, options | subs_options::no_pattern);
606 // avoid infinite recursion when re-substituting the wildcards
613 /** Substitute a set of objects by arbitrary expressions. The ex returned
614 * will already be evaluated. */
615 ex basic::subs(const exmap & m, unsigned options) const
620 // Substitute in subexpressions
621 for (size_t i=0; i<num; i++) {
622 const ex & orig_op = op(i);
623 const ex & subsed_op = orig_op.subs(m, options);
624 if (!are_ex_trivially_equal(orig_op, subsed_op)) {
626 // Something changed, clone the object
627 basic *copy = duplicate();
628 copy->setflag(status_flags::dynallocated);
629 copy->clearflag(status_flags::hash_calculated | status_flags::expanded);
631 // Substitute the changed operand
632 copy->let_op(i++) = subsed_op;
634 // Substitute the other operands
636 copy->let_op(i) = op(i).subs(m, options);
638 // Perform substitutions on the new object as a whole
639 return copy->subs_one_level(m, options);
644 // Nothing changed or no subexpressions
645 return subs_one_level(m, options);
648 /** Default interface of nth derivative ex::diff(s, n). It should be called
649 * instead of ::derivative(s) for first derivatives and for nth derivatives it
650 * just recurses down.
652 * @param s symbol to differentiate in
653 * @param nth order of differentiation
655 ex basic::diff(const symbol & s, unsigned nth) const
657 // trivial: zeroth derivative
661 // evaluate unevaluated *this before differentiating
662 if (!(flags & status_flags::evaluated))
663 return ex(*this).diff(s, nth);
665 ex ndiff = this->derivative(s);
666 while (!ndiff.is_zero() && // stop differentiating zeros
668 ndiff = ndiff.diff(s);
674 /** Return a vector containing the free indices of an expression. */
675 exvector basic::get_free_indices() const
677 return exvector(); // return an empty exvector
680 ex basic::conjugate() const
685 ex basic::real_part() const
687 return real_part_function(*this).hold();
690 ex basic::imag_part() const
692 return imag_part_function(*this).hold();
695 ex basic::eval_ncmul(const exvector & v) const
697 return hold_ncmul(v);
702 /** Function object to be applied by basic::derivative(). */
703 struct derivative_map_function : public map_function {
705 derivative_map_function(const symbol &sym) : s(sym) {}
706 ex operator()(const ex & e) { return diff(e, s); }
709 /** Default implementation of ex::diff(). It maps the operation on the
710 * operands (or returns 0 when the object has no operands).
713 ex basic::derivative(const symbol & s) const
718 derivative_map_function map_derivative(s);
719 return map(map_derivative);
723 /** Returns order relation between two objects of same type. This needs to be
724 * implemented by each class. It may never return anything else than 0,
725 * signalling equality, or +1 and -1 signalling inequality and determining
726 * the canonical ordering. (Perl hackers will wonder why C++ doesn't feature
727 * the spaceship operator <=> for denoting just this.) */
728 int basic::compare_same_type(const basic & other) const
730 return compare_pointers(this, &other);
733 /** Returns true if two objects of same type are equal. Normally needs
734 * not be reimplemented as long as it wasn't overwritten by some parent
735 * class, since it just calls compare_same_type(). The reason why this
736 * function exists is that sometimes it is easier to determine equality
737 * than an order relation and then it can be overridden. */
738 bool basic::is_equal_same_type(const basic & other) const
740 return compare_same_type(other)==0;
743 /** Returns true if the attributes of two objects are similar enough for
744 * a match. This function must not match subexpressions (this is already
745 * done by basic::match()). Only attributes not accessible by op() should
746 * be compared. This is also the reason why this function doesn't take the
747 * wildcard replacement list from match() as an argument: only subexpressions
748 * are subject to wildcard matches. Also, this function only needs to be
749 * implemented for container classes because is_equal_same_type() is
750 * automatically used instead of match_same_type() if nops() == 0.
752 * @see basic::match */
753 bool basic::match_same_type(const basic & other) const
755 // The default is to only consider subexpressions, but not any other
760 unsigned basic::return_type() const
762 return return_types::commutative;
765 return_type_t basic::return_type_tinfo() const
768 rt.tinfo = &typeid(*this);
773 /** Compute the hash value of an object and if it makes sense to store it in
774 * the objects status_flags, do so. The method inherited from class basic
775 * computes a hash value based on the type and hash values of possible
776 * members. For this reason it is well suited for container classes but
777 * atomic classes should override this implementation because otherwise they
778 * would all end up with the same hashvalue. */
779 unsigned basic::calchash() const
781 unsigned v = make_hash_seed(typeid(*this));
782 for (size_t i=0; i<nops(); i++) {
784 v ^= this->op(i).gethash();
787 // store calculated hash value only if object is already evaluated
788 if (flags & status_flags::evaluated) {
789 setflag(status_flags::hash_calculated);
796 /** Function object to be applied by basic::expand(). */
797 struct expand_map_function : public map_function {
799 expand_map_function(unsigned o) : options(o) {}
800 ex operator()(const ex & e) { return e.expand(options); }
803 /** Expand expression, i.e. multiply it out and return the result as a new
805 ex basic::expand(unsigned options) const
808 return (options == 0) ? setflag(status_flags::expanded) : *this;
810 expand_map_function map_expand(options);
811 return ex_to<basic>(map(map_expand)).setflag(options == 0 ? status_flags::expanded : 0);
817 // non-virtual functions in this class
822 /** Compare objects syntactically to establish canonical ordering.
823 * All compare functions return: -1 for *this less than other, 0 equal,
825 int basic::compare(const basic & other) const
827 #ifdef GINAC_COMPARE_STATISTICS
828 compare_statistics.total_basic_compares++;
830 const unsigned hash_this = gethash();
831 const unsigned hash_other = other.gethash();
832 if (hash_this<hash_other) return -1;
833 if (hash_this>hash_other) return 1;
834 #ifdef GINAC_COMPARE_STATISTICS
835 compare_statistics.compare_same_hashvalue++;
838 const std::type_info& typeid_this = typeid(*this);
839 const std::type_info& typeid_other = typeid(other);
840 if (typeid_this == typeid_other) {
841 // int cmpval = compare_same_type(other);
843 // std::cout << "hash collision, same type: "
844 // << *this << " and " << other << std::endl;
845 // this->print(print_tree(std::cout));
846 // std::cout << " and ";
847 // other.print(print_tree(std::cout));
848 // std::cout << std::endl;
851 #ifdef GINAC_COMPARE_STATISTICS
852 compare_statistics.compare_same_type++;
854 return compare_same_type(other);
856 // std::cout << "hash collision, different types: "
857 // << *this << " and " << other << std::endl;
858 // this->print(print_tree(std::cout));
859 // std::cout << " and ";
860 // other.print(print_tree(std::cout));
861 // std::cout << std::endl;
862 return (typeid_this.before(typeid_other) ? -1 : 1);
866 /** Test for syntactic equality.
867 * This is only a quick test, meaning objects should be in the same domain.
868 * You might have to .expand(), .normal() objects first, depending on the
869 * domain of your computation, to get a more reliable answer.
871 * @see is_equal_same_type */
872 bool basic::is_equal(const basic & other) const
874 #ifdef GINAC_COMPARE_STATISTICS
875 compare_statistics.total_basic_is_equals++;
877 if (this->gethash()!=other.gethash())
879 #ifdef GINAC_COMPARE_STATISTICS
880 compare_statistics.is_equal_same_hashvalue++;
882 if (typeid(*this) != typeid(other))
885 #ifdef GINAC_COMPARE_STATISTICS
886 compare_statistics.is_equal_same_type++;
888 return is_equal_same_type(other);
893 /** Stop further evaluation.
895 * @see basic::eval */
896 const basic & basic::hold() const
898 return setflag(status_flags::evaluated);
901 /** Ensure the object may be modified without hurting others, throws if this
902 * is not the case. */
903 void basic::ensure_if_modifiable() const
905 if (get_refcount() > 1)
906 throw(std::runtime_error("cannot modify multiply referenced object"));
907 clearflag(status_flags::hash_calculated | status_flags::evaluated);
914 int max_recursion_level = 1024;
917 #ifdef GINAC_COMPARE_STATISTICS
918 compare_statistics_t::~compare_statistics_t()
920 std::clog << "ex::compare() called " << total_compares << " times" << std::endl;
921 std::clog << "nontrivial compares: " << nontrivial_compares << " times" << std::endl;
922 std::clog << "basic::compare() called " << total_basic_compares << " times" << std::endl;
923 std::clog << "same hashvalue in compare(): " << compare_same_hashvalue << " times" << std::endl;
924 std::clog << "compare_same_type() called " << compare_same_type << " times" << std::endl;
925 std::clog << std::endl;
926 std::clog << "ex::is_equal() called " << total_is_equals << " times" << std::endl;
927 std::clog << "nontrivial is_equals: " << nontrivial_is_equals << " times" << std::endl;
928 std::clog << "basic::is_equal() called " << total_basic_is_equals << " times" << std::endl;
929 std::clog << "same hashvalue in is_equal(): " << is_equal_same_hashvalue << " times" << std::endl;
930 std::clog << "is_equal_same_type() called " << is_equal_same_type << " times" << std::endl;
931 std::clog << std::endl;
932 std::clog << "basic::gethash() called " << total_gethash << " times" << std::endl;
933 std::clog << "used cached hashvalue " << gethash_cached << " times" << std::endl;
936 compare_statistics_t compare_statistics;