3 * Implementation of GiNaC's ABC. */
6 * GiNaC Copyright (C) 1999-2003 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., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
25 #ifdef DO_GINAC_ASSERT
36 #include "relational.h"
37 #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), refcount(0)
60 GINAC_ASSERT(typeid(*this) == typeid(other));
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 (tinfo_key != other.tinfo_key) {
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.
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 basic::basic(const archive_node &n, lst &sym_lst) : flags(0), refcount(0)
98 // Reconstruct tinfo_key from class name
99 std::string class_name;
100 if (n.find_string("class", class_name))
101 tinfo_key = find_tinfo_key(class_name);
103 throw (std::runtime_error("archive node contains no class name"));
106 /** Unarchive the object. */
107 DEFAULT_UNARCHIVE(basic)
109 /** Archive the object. */
110 void basic::archive(archive_node &n) const
112 n.add_string("class", class_name());
116 // new virtual functions which can be overridden by derived classes
121 /** Output to stream. This performs double dispatch on the dynamic type of
122 * *this and the dynamic type of the supplied print context.
123 * @param c print context object that describes the output formatting
124 * @param level value that is used to identify the precedence or indentation
125 * level for placing parentheses and formatting */
126 void basic::print(const print_context & c, unsigned level) const
128 print_dispatch(get_class_info(), c, level);
131 /** Like print(), but dispatch to the specified class. Can be used by
132 * implementations of print methods to dispatch to the method of the
135 * @see basic::print */
136 void basic::print_dispatch(const registered_class_info & ri, const print_context & c, unsigned level) const
138 // Double dispatch on object type and print_context type
139 const registered_class_info * reg_info = &ri;
140 const print_context_class_info * pc_info = &c.get_class_info();
143 const std::vector<print_functor> & pdt = reg_info->options.get_print_dispatch_table();
146 unsigned id = pc_info->options.get_id();
147 if (id >= pdt.size() || !(pdt[id].is_valid())) {
149 // Method not found, try parent print_context class
150 const print_context_class_info * parent_pc_info = pc_info->get_parent();
151 if (parent_pc_info) {
152 pc_info = parent_pc_info;
156 // Method still not found, try parent class
157 const registered_class_info * parent_reg_info = reg_info->get_parent();
158 if (parent_reg_info) {
159 reg_info = parent_reg_info;
160 pc_info = &c.get_class_info();
164 // Method still not found. This shouldn't happen because basic (the
165 // base class of the algebraic hierarchy) registers a method for
166 // print_context (the base class of the print context hierarchy),
167 // so if we end up here, there's something wrong with the class
169 throw (std::runtime_error(std::string("basic::print(): method for ") + class_name() + "/" + c.class_name() + " not found"));
174 pdt[id](*this, c, level);
178 /** Default output to stream. */
179 void basic::do_print(const print_context & c, unsigned level) const
181 c.s << "[" << class_name() << " object]";
184 /** Tree output to stream. */
185 void basic::do_print_tree(const print_tree & c, unsigned level) const
187 c.s << std::string(level, ' ') << class_name()
188 << std::hex << ", hash=0x" << hashvalue << ", flags=0x" << flags << std::dec;
190 c.s << ", nops=" << nops();
192 for (size_t i=0; i<nops(); ++i)
193 op(i).print(c, level + c.delta_indent);
196 /** Python parsable output to stream. */
197 void basic::do_print_python_repr(const print_python_repr & c, unsigned level) const
199 c.s << class_name() << "()";
202 /** Little wrapper around print to be called within a debugger.
203 * This is needed because you cannot call foo.print(cout) from within the
204 * debugger because it might not know what cout is. This method can be
205 * invoked with no argument and it will simply print to stdout.
208 * @see basic::dbgprinttree */
209 void basic::dbgprint() const
211 this->print(std::cerr);
212 std::cerr << std::endl;
215 /** Little wrapper around printtree to be called within a debugger.
217 * @see basic::dbgprint */
218 void basic::dbgprinttree() const
220 this->print(print_tree(std::cerr));
223 /** Return relative operator precedence (for parenthezing output). */
224 unsigned basic::precedence() const
229 /** Information about the object.
231 * @see class info_flags */
232 bool basic::info(unsigned inf) const
234 // all possible properties are false for basic objects
238 /** Number of operands/members. */
239 size_t basic::nops() const
241 // iterating from 0 to nops() on atomic objects should be an empty loop,
242 // and accessing their elements is a range error. Container objects should
247 /** Return operand/member at position i. */
248 ex basic::op(size_t i) const
250 throw(std::range_error(std::string("basic::op(): ") + class_name() + std::string(" has no operands")));
253 /** Return modifyable operand/member at position i. */
254 ex & basic::let_op(size_t i)
256 ensure_if_modifiable();
257 throw(std::range_error(std::string("basic::let_op(): ") + class_name() + std::string(" has no operands")));
260 ex basic::operator[](const ex & index) const
262 if (is_exactly_a<numeric>(index))
263 return op(static_cast<size_t>(ex_to<numeric>(index).to_int()));
265 throw(std::invalid_argument(std::string("non-numeric indices not supported by ") + class_name()));
268 ex basic::operator[](size_t i) const
273 ex & basic::operator[](const ex & index)
275 if (is_exactly_a<numeric>(index))
276 return let_op(ex_to<numeric>(index).to_int());
278 throw(std::invalid_argument(std::string("non-numeric indices not supported by ") + class_name()));
281 ex & basic::operator[](size_t i)
286 /** Test for occurrence of a pattern. An object 'has' a pattern if it matches
287 * the pattern itself or one of the children 'has' it. As a consequence
288 * (according to the definition of children) given e=x+y+z, e.has(x) is true
289 * but e.has(x+y) is false. */
290 bool basic::has(const ex & pattern) const
293 if (match(pattern, repl_lst))
295 for (size_t i=0; i<nops(); i++)
296 if (op(i).has(pattern))
302 /** Construct new expression by applying the specified function to all
303 * sub-expressions (one level only, not recursively). */
304 ex basic::map(map_function & f) const
310 basic *copy = duplicate();
311 copy->setflag(status_flags::dynallocated);
312 copy->clearflag(status_flags::hash_calculated | status_flags::expanded);
313 for (size_t i=0; i<num; i++)
314 copy->let_op(i) = f(copy->op(i));
318 /** Return degree of highest power in object s. */
319 int basic::degree(const ex & s) const
321 return is_equal(ex_to<basic>(s)) ? 1 : 0;
324 /** Return degree of lowest power in object s. */
325 int basic::ldegree(const ex & s) const
327 return is_equal(ex_to<basic>(s)) ? 1 : 0;
330 /** Return coefficient of degree n in object s. */
331 ex basic::coeff(const ex & s, int n) const
333 if (is_equal(ex_to<basic>(s)))
334 return n==1 ? _ex1 : _ex0;
336 return n==0 ? *this : _ex0;
339 /** Sort expanded expression in terms of powers of some object(s).
340 * @param s object(s) to sort in
341 * @param distributed recursive or distributed form (only used when s is a list) */
342 ex basic::collect(const ex & s, bool distributed) const
347 // List of objects specified
351 return collect(s.op(0));
353 else if (distributed) {
355 // Get lower/upper degree of all symbols in list
356 size_t num = s.nops();
360 int cnt; // current degree, 'counter'
361 ex coeff; // coefficient for degree 'cnt'
363 sym_info *si = new sym_info[num];
365 for (size_t i=0; i<num; i++) {
367 si[i].ldeg = si[i].cnt = this->ldegree(si[i].sym);
368 si[i].deg = this->degree(si[i].sym);
369 c = si[i].coeff = c.coeff(si[i].sym, si[i].cnt);
374 // Calculate coeff*x1^c1*...*xn^cn
376 for (size_t i=0; i<num; i++) {
378 y *= power(si[i].sym, cnt);
380 x += y * si[num - 1].coeff;
382 // Increment counters
386 if (si[n].cnt <= si[n].deg) {
387 // Update coefficients
393 for (size_t i=n; i<num; i++)
394 c = si[i].coeff = c.coeff(si[i].sym, si[i].cnt);
399 si[n].cnt = si[n].ldeg;
410 size_t n = s.nops() - 1;
421 // Only one object specified
422 for (int n=this->ldegree(s); n<=this->degree(s); ++n)
423 x += this->coeff(s,n)*power(s,n);
426 // correct for lost fractional arguments and return
427 return x + (*this - x).expand();
430 /** Perform automatic non-interruptive term rewriting rules. */
431 ex basic::eval(int level) const
433 // There is nothing to do for basic objects:
437 /** Function object to be applied by basic::evalf(). */
438 struct evalf_map_function : public map_function {
440 evalf_map_function(int l) : level(l) {}
441 ex operator()(const ex & e) { return evalf(e, level); }
444 /** Evaluate object numerically. */
445 ex basic::evalf(int level) const
452 else if (level == -max_recursion_level)
453 throw(std::runtime_error("max recursion level reached"));
455 evalf_map_function map_evalf(level - 1);
456 return map(map_evalf);
461 /** Function object to be applied by basic::evalm(). */
462 struct evalm_map_function : public map_function {
463 ex operator()(const ex & e) { return evalm(e); }
466 /** Evaluate sums, products and integer powers of matrices. */
467 ex basic::evalm() const
472 return map(map_evalm);
475 /** Perform automatic symbolic evaluations on indexed expression that
476 * contains this object as the base expression. */
477 ex basic::eval_indexed(const basic & i) const
478 // this function can't take a "const ex & i" because that would result
479 // in an infinite eval() loop
481 // There is nothing to do for basic objects
485 /** Add two indexed expressions. They are guaranteed to be of class indexed
486 * (or a subclass) and their indices are compatible. This function is used
487 * internally by simplify_indexed().
489 * @param self First indexed expression; it's base object is *this
490 * @param other Second indexed expression
491 * @return sum of self and other
492 * @see ex::simplify_indexed() */
493 ex basic::add_indexed(const ex & self, const ex & other) const
498 /** Multiply an indexed expression with a scalar. This function is used
499 * internally by simplify_indexed().
501 * @param self Indexed expression; it's base object is *this
502 * @param other Numeric value
503 * @return product of self and other
504 * @see ex::simplify_indexed() */
505 ex basic::scalar_mul_indexed(const ex & self, const numeric & other) const
510 /** Try to contract two indexed expressions that appear in the same product.
511 * If a contraction exists, the function overwrites one or both of the
512 * expressions and returns true. Otherwise it returns false. It is
513 * guaranteed that both expressions are of class indexed (or a subclass)
514 * and that at least one dummy index has been found. This functions is
515 * used internally by simplify_indexed().
517 * @param self Pointer to first indexed expression; it's base object is *this
518 * @param other Pointer to second indexed expression
519 * @param v The complete vector of factors
520 * @return true if the contraction was successful, false otherwise
521 * @see ex::simplify_indexed() */
522 bool basic::contract_with(exvector::iterator self, exvector::iterator other, exvector & v) const
528 /** Check whether the expression matches a given pattern. For every wildcard
529 * object in the pattern, an expression of the form "wildcard == matching_expression"
530 * is added to repl_lst. */
531 bool basic::match(const ex & pattern, lst & repl_lst) const
534 Sweet sweet shapes, sweet sweet shapes,
535 That's the key thing, right right.
536 Feed feed face, feed feed shapes,
537 But who is the king tonight?
538 Who is the king tonight?
539 Pattern is the thing, the key thing-a-ling,
540 But who is the king of Pattern?
541 But who is the king, the king thing-a-ling,
542 Who is the king of Pattern?
543 Bog is the king, the king thing-a-ling,
544 Bog is the king of Pattern.
545 Ba bu-bu-bu-bu bu-bu-bu-bu-bu-bu bu-bu
546 Bog is the king of Pattern.
549 if (is_exactly_a<wildcard>(pattern)) {
551 // Wildcard matches anything, but check whether we already have found
552 // a match for that wildcard first (if so, it the earlier match must
553 // be the same expression)
554 for (lst::const_iterator it = repl_lst.begin(); it != repl_lst.end(); ++it) {
555 if (it->op(0).is_equal(pattern))
556 return is_equal(ex_to<basic>(it->op(1)));
558 repl_lst.append(pattern == *this);
563 // Expression must be of the same type as the pattern
564 if (tinfo() != ex_to<basic>(pattern).tinfo())
567 // Number of subexpressions must match
568 if (nops() != pattern.nops())
571 // No subexpressions? Then just compare the objects (there can't be
572 // wildcards in the pattern)
574 return is_equal_same_type(ex_to<basic>(pattern));
576 // Check whether attributes that are not subexpressions match
577 if (!match_same_type(ex_to<basic>(pattern)))
580 // Otherwise the subexpressions must match one-to-one
581 for (size_t i=0; i<nops(); i++)
582 if (!op(i).match(pattern.op(i), repl_lst))
585 // Looks similar enough, match found
590 /** Helper function for subs(). Does not recurse into subexpressions. */
591 ex basic::subs_one_level(const exmap & m, unsigned options) const
593 exmap::const_iterator it;
595 if (options & subs_options::no_pattern) {
600 for (it = m.begin(); it != m.end(); ++it) {
602 if (match(ex_to<basic>(it->first), repl_lst))
603 return it->second.subs(repl_lst, options | subs_options::no_pattern); // avoid infinite recursion when re-substituting the wildcards
610 /** Substitute a set of objects by arbitrary expressions. The ex returned
611 * will already be evaluated. */
612 ex basic::subs(const exmap & m, unsigned options) const
617 // Substitute in subexpressions
618 for (size_t i=0; i<num; i++) {
619 const ex & orig_op = op(i);
620 const ex & subsed_op = orig_op.subs(m, options);
621 if (!are_ex_trivially_equal(orig_op, subsed_op)) {
623 // Something changed, clone the object
624 basic *copy = duplicate();
625 copy->setflag(status_flags::dynallocated);
626 copy->clearflag(status_flags::hash_calculated | status_flags::expanded);
628 // Substitute the changed operand
629 copy->let_op(i++) = subsed_op;
631 // Substitute the other operands
633 copy->let_op(i) = op(i).subs(m, options);
635 // Perform substitutions on the new object as a whole
636 return copy->subs_one_level(m, options);
641 // Nothing changed or no subexpressions
642 return subs_one_level(m, options);
645 /** Default interface of nth derivative ex::diff(s, n). It should be called
646 * instead of ::derivative(s) for first derivatives and for nth derivatives it
647 * just recurses down.
649 * @param s symbol to differentiate in
650 * @param nth order of differentiation
652 ex basic::diff(const symbol & s, unsigned nth) const
654 // trivial: zeroth derivative
658 // evaluate unevaluated *this before differentiating
659 if (!(flags & status_flags::evaluated))
660 return ex(*this).diff(s, nth);
662 ex ndiff = this->derivative(s);
663 while (!ndiff.is_zero() && // stop differentiating zeros
665 ndiff = ndiff.diff(s);
671 /** Return a vector containing the free indices of an expression. */
672 exvector basic::get_free_indices() const
674 return exvector(); // return an empty exvector
677 ex basic::eval_ncmul(const exvector & v) const
679 return hold_ncmul(v);
684 /** Function object to be applied by basic::derivative(). */
685 struct derivative_map_function : public map_function {
687 derivative_map_function(const symbol &sym) : s(sym) {}
688 ex operator()(const ex & e) { return diff(e, s); }
691 /** Default implementation of ex::diff(). It maps the operation on the
692 * operands (or returns 0 when the object has no operands).
695 ex basic::derivative(const symbol & s) const
700 derivative_map_function map_derivative(s);
701 return map(map_derivative);
705 /** Returns order relation between two objects of same type. This needs to be
706 * implemented by each class. It may never return anything else than 0,
707 * signalling equality, or +1 and -1 signalling inequality and determining
708 * the canonical ordering. (Perl hackers will wonder why C++ doesn't feature
709 * the spaceship operator <=> for denoting just this.) */
710 int basic::compare_same_type(const basic & other) const
712 return compare_pointers(this, &other);
715 /** Returns true if two objects of same type are equal. Normally needs
716 * not be reimplemented as long as it wasn't overwritten by some parent
717 * class, since it just calls compare_same_type(). The reason why this
718 * function exists is that sometimes it is easier to determine equality
719 * than an order relation and then it can be overridden. */
720 bool basic::is_equal_same_type(const basic & other) const
722 return compare_same_type(other)==0;
725 /** Returns true if the attributes of two objects are similar enough for
726 * a match. This function must not match subexpressions (this is already
727 * done by basic::match()). Only attributes not accessible by op() should
728 * be compared. This is also the reason why this function doesn't take the
729 * wildcard replacement list from match() as an argument: only subexpressions
730 * are subject to wildcard matches. Also, this function only needs to be
731 * implemented for container classes because is_equal_same_type() is
732 * automatically used instead of match_same_type() if nops() == 0.
734 * @see basic::match */
735 bool basic::match_same_type(const basic & other) const
737 // The default is to only consider subexpressions, but not any other
742 unsigned basic::return_type() const
744 return return_types::commutative;
747 unsigned basic::return_type_tinfo() const
752 /** Compute the hash value of an object and if it makes sense to store it in
753 * the objects status_flags, do so. The method inherited from class basic
754 * computes a hash value based on the type and hash values of possible
755 * members. For this reason it is well suited for container classes but
756 * atomic classes should override this implementation because otherwise they
757 * would all end up with the same hashvalue. */
758 unsigned basic::calchash() const
760 unsigned v = golden_ratio_hash(tinfo());
761 for (size_t i=0; i<nops(); i++) {
763 v ^= this->op(i).gethash();
766 // store calculated hash value only if object is already evaluated
767 if (flags & status_flags::evaluated) {
768 setflag(status_flags::hash_calculated);
775 /** Function object to be applied by basic::expand(). */
776 struct expand_map_function : public map_function {
778 expand_map_function(unsigned o) : options(o) {}
779 ex operator()(const ex & e) { return expand(e, options); }
782 /** Expand expression, i.e. multiply it out and return the result as a new
784 ex basic::expand(unsigned options) const
787 return (options == 0) ? setflag(status_flags::expanded) : *this;
789 expand_map_function map_expand(options);
790 return ex_to<basic>(map(map_expand)).setflag(options == 0 ? status_flags::expanded : 0);
796 // non-virtual functions in this class
801 /** Compare objects syntactically to establish canonical ordering.
802 * All compare functions return: -1 for *this less than other, 0 equal,
804 int basic::compare(const basic & other) const
806 const unsigned hash_this = gethash();
807 const unsigned hash_other = other.gethash();
808 if (hash_this<hash_other) return -1;
809 if (hash_this>hash_other) return 1;
811 const unsigned typeid_this = tinfo();
812 const unsigned typeid_other = other.tinfo();
813 if (typeid_this==typeid_other) {
814 GINAC_ASSERT(typeid(*this)==typeid(other));
815 // int cmpval = compare_same_type(other);
817 // std::cout << "hash collision, same type: "
818 // << *this << " and " << other << std::endl;
819 // this->print(print_tree(std::cout));
820 // std::cout << " and ";
821 // other.print(print_tree(std::cout));
822 // std::cout << std::endl;
825 return compare_same_type(other);
827 // std::cout << "hash collision, different types: "
828 // << *this << " and " << other << std::endl;
829 // this->print(print_tree(std::cout));
830 // std::cout << " and ";
831 // other.print(print_tree(std::cout));
832 // std::cout << std::endl;
833 return (typeid_this<typeid_other ? -1 : 1);
837 /** Test for syntactic equality.
838 * This is only a quick test, meaning objects should be in the same domain.
839 * You might have to .expand(), .normal() objects first, depending on the
840 * domain of your computation, to get a more reliable answer.
842 * @see is_equal_same_type */
843 bool basic::is_equal(const basic & other) const
845 if (this->gethash()!=other.gethash())
847 if (this->tinfo()!=other.tinfo())
850 GINAC_ASSERT(typeid(*this)==typeid(other));
852 return is_equal_same_type(other);
857 /** Stop further evaluation.
859 * @see basic::eval */
860 const basic & basic::hold() const
862 return setflag(status_flags::evaluated);
865 /** Ensure the object may be modified without hurting others, throws if this
866 * is not the case. */
867 void basic::ensure_if_modifiable() const
870 throw(std::runtime_error("cannot modify multiply referenced object"));
871 clearflag(status_flags::hash_calculated | status_flags::evaluated);
878 int max_recursion_level = 1024;