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
189 << ", nops=" << nops()
191 for (size_t i=0; i<nops(); ++i)
192 op(i).print(c, level + c.delta_indent);
195 /** Python parsable output to stream. */
196 void basic::do_print_python_repr(const print_python_repr & c, unsigned level) const
198 c.s << class_name() << "()";
201 /** Little wrapper around print to be called within a debugger.
202 * This is needed because you cannot call foo.print(cout) from within the
203 * debugger because it might not know what cout is. This method can be
204 * invoked with no argument and it will simply print to stdout.
207 * @see basic::dbgprinttree */
208 void basic::dbgprint() const
210 this->print(std::cerr);
211 std::cerr << std::endl;
214 /** Little wrapper around printtree to be called within a debugger.
216 * @see basic::dbgprint */
217 void basic::dbgprinttree() const
219 this->print(print_tree(std::cerr));
222 /** Return relative operator precedence (for parenthezing output). */
223 unsigned basic::precedence() const
228 /** Information about the object.
230 * @see class info_flags */
231 bool basic::info(unsigned inf) const
233 // all possible properties are false for basic objects
237 /** Number of operands/members. */
238 size_t basic::nops() const
240 // iterating from 0 to nops() on atomic objects should be an empty loop,
241 // and accessing their elements is a range error. Container objects should
246 /** Return operand/member at position i. */
247 ex basic::op(size_t i) const
249 throw(std::range_error(std::string("basic::op(): ") + class_name() + std::string(" has no operands")));
252 /** Return modifyable operand/member at position i. */
253 ex & basic::let_op(size_t i)
255 ensure_if_modifiable();
256 throw(std::range_error(std::string("basic::let_op(): ") + class_name() + std::string(" has no operands")));
259 ex basic::operator[](const ex & index) const
261 if (is_exactly_a<numeric>(index))
262 return op(static_cast<size_t>(ex_to<numeric>(index).to_int()));
264 throw(std::invalid_argument(std::string("non-numeric indices not supported by ") + class_name()));
267 ex basic::operator[](size_t i) const
272 ex & basic::operator[](const ex & index)
274 if (is_exactly_a<numeric>(index))
275 return let_op(ex_to<numeric>(index).to_int());
277 throw(std::invalid_argument(std::string("non-numeric indices not supported by ") + class_name()));
280 ex & basic::operator[](size_t i)
285 /** Test for occurrence of a pattern. An object 'has' a pattern if it matches
286 * the pattern itself or one of the children 'has' it. As a consequence
287 * (according to the definition of children) given e=x+y+z, e.has(x) is true
288 * but e.has(x+y) is false. */
289 bool basic::has(const ex & pattern) const
292 if (match(pattern, repl_lst))
294 for (size_t i=0; i<nops(); i++)
295 if (op(i).has(pattern))
301 /** Construct new expression by applying the specified function to all
302 * sub-expressions (one level only, not recursively). */
303 ex basic::map(map_function & f) const
309 basic *copy = duplicate();
310 copy->setflag(status_flags::dynallocated);
311 copy->clearflag(status_flags::hash_calculated | status_flags::expanded);
312 for (size_t i=0; i<num; i++)
313 copy->let_op(i) = f(copy->op(i));
317 /** Return degree of highest power in object s. */
318 int basic::degree(const ex & s) const
320 return is_equal(ex_to<basic>(s)) ? 1 : 0;
323 /** Return degree of lowest power in object s. */
324 int basic::ldegree(const ex & s) const
326 return is_equal(ex_to<basic>(s)) ? 1 : 0;
329 /** Return coefficient of degree n in object s. */
330 ex basic::coeff(const ex & s, int n) const
332 if (is_equal(ex_to<basic>(s)))
333 return n==1 ? _ex1 : _ex0;
335 return n==0 ? *this : _ex0;
338 /** Sort expanded expression in terms of powers of some object(s).
339 * @param s object(s) to sort in
340 * @param distributed recursive or distributed form (only used when s is a list) */
341 ex basic::collect(const ex & s, bool distributed) const
346 // List of objects specified
350 return collect(s.op(0));
352 else if (distributed) {
354 // Get lower/upper degree of all symbols in list
355 size_t num = s.nops();
359 int cnt; // current degree, 'counter'
360 ex coeff; // coefficient for degree 'cnt'
362 sym_info *si = new sym_info[num];
364 for (size_t i=0; i<num; i++) {
366 si[i].ldeg = si[i].cnt = this->ldegree(si[i].sym);
367 si[i].deg = this->degree(si[i].sym);
368 c = si[i].coeff = c.coeff(si[i].sym, si[i].cnt);
373 // Calculate coeff*x1^c1*...*xn^cn
375 for (size_t i=0; i<num; i++) {
377 y *= power(si[i].sym, cnt);
379 x += y * si[num - 1].coeff;
381 // Increment counters
385 if (si[n].cnt <= si[n].deg) {
386 // Update coefficients
392 for (size_t i=n; i<num; i++)
393 c = si[i].coeff = c.coeff(si[i].sym, si[i].cnt);
398 si[n].cnt = si[n].ldeg;
409 size_t n = s.nops() - 1;
420 // Only one object specified
421 for (int n=this->ldegree(s); n<=this->degree(s); ++n)
422 x += this->coeff(s,n)*power(s,n);
425 // correct for lost fractional arguments and return
426 return x + (*this - x).expand();
429 /** Perform automatic non-interruptive term rewriting rules. */
430 ex basic::eval(int level) const
432 // There is nothing to do for basic objects:
436 /** Function object to be applied by basic::evalf(). */
437 struct evalf_map_function : public map_function {
439 evalf_map_function(int l) : level(l) {}
440 ex operator()(const ex & e) { return evalf(e, level); }
443 /** Evaluate object numerically. */
444 ex basic::evalf(int level) const
451 else if (level == -max_recursion_level)
452 throw(std::runtime_error("max recursion level reached"));
454 evalf_map_function map_evalf(level - 1);
455 return map(map_evalf);
460 /** Function object to be applied by basic::evalm(). */
461 struct evalm_map_function : public map_function {
462 ex operator()(const ex & e) { return evalm(e); }
465 /** Evaluate sums, products and integer powers of matrices. */
466 ex basic::evalm() const
471 return map(map_evalm);
474 /** Perform automatic symbolic evaluations on indexed expression that
475 * contains this object as the base expression. */
476 ex basic::eval_indexed(const basic & i) const
477 // this function can't take a "const ex & i" because that would result
478 // in an infinite eval() loop
480 // There is nothing to do for basic objects
484 /** Add two indexed expressions. They are guaranteed to be of class indexed
485 * (or a subclass) and their indices are compatible. This function is used
486 * internally by simplify_indexed().
488 * @param self First indexed expression; it's base object is *this
489 * @param other Second indexed expression
490 * @return sum of self and other
491 * @see ex::simplify_indexed() */
492 ex basic::add_indexed(const ex & self, const ex & other) const
497 /** Multiply an indexed expression with a scalar. This function is used
498 * internally by simplify_indexed().
500 * @param self Indexed expression; it's base object is *this
501 * @param other Numeric value
502 * @return product of self and other
503 * @see ex::simplify_indexed() */
504 ex basic::scalar_mul_indexed(const ex & self, const numeric & other) const
509 /** Try to contract two indexed expressions that appear in the same product.
510 * If a contraction exists, the function overwrites one or both of the
511 * expressions and returns true. Otherwise it returns false. It is
512 * guaranteed that both expressions are of class indexed (or a subclass)
513 * and that at least one dummy index has been found. This functions is
514 * used internally by simplify_indexed().
516 * @param self Pointer to first indexed expression; it's base object is *this
517 * @param other Pointer to second indexed expression
518 * @param v The complete vector of factors
519 * @return true if the contraction was successful, false otherwise
520 * @see ex::simplify_indexed() */
521 bool basic::contract_with(exvector::iterator self, exvector::iterator other, exvector & v) const
527 /** Check whether the expression matches a given pattern. For every wildcard
528 * object in the pattern, an expression of the form "wildcard == matching_expression"
529 * is added to repl_lst. */
530 bool basic::match(const ex & pattern, lst & repl_lst) const
533 Sweet sweet shapes, sweet sweet shapes,
534 That's the key thing, right right.
535 Feed feed face, feed feed shapes,
536 But who is the king tonight?
537 Who is the king tonight?
538 Pattern is the thing, the key thing-a-ling,
539 But who is the king of Pattern?
540 But who is the king, the king thing-a-ling,
541 Who is the king of Pattern?
542 Bog is the king, the king thing-a-ling,
543 Bog is the king of Pattern.
544 Ba bu-bu-bu-bu bu-bu-bu-bu-bu-bu bu-bu
545 Bog is the king of Pattern.
548 if (is_exactly_a<wildcard>(pattern)) {
550 // Wildcard matches anything, but check whether we already have found
551 // a match for that wildcard first (if so, it the earlier match must
552 // be the same expression)
553 for (lst::const_iterator it = repl_lst.begin(); it != repl_lst.end(); ++it) {
554 if (it->op(0).is_equal(pattern))
555 return is_equal(ex_to<basic>(it->op(1)));
557 repl_lst.append(pattern == *this);
562 // Expression must be of the same type as the pattern
563 if (tinfo() != ex_to<basic>(pattern).tinfo())
566 // Number of subexpressions must match
567 if (nops() != pattern.nops())
570 // No subexpressions? Then just compare the objects (there can't be
571 // wildcards in the pattern)
573 return is_equal_same_type(ex_to<basic>(pattern));
575 // Check whether attributes that are not subexpressions match
576 if (!match_same_type(ex_to<basic>(pattern)))
579 // Otherwise the subexpressions must match one-to-one
580 for (size_t i=0; i<nops(); i++)
581 if (!op(i).match(pattern.op(i), repl_lst))
584 // Looks similar enough, match found
589 /** Helper function for subs(). Does not recurse into subexpressions. */
590 ex basic::subs_one_level(const exmap & m, unsigned options) const
592 exmap::const_iterator it;
594 if (options & subs_options::no_pattern) {
599 for (it = m.begin(); it != m.end(); ++it) {
601 if (match(ex_to<basic>(it->first), repl_lst))
602 return it->second.subs(repl_lst, options | subs_options::no_pattern); // avoid infinite recursion when re-substituting the wildcards
609 /** Substitute a set of objects by arbitrary expressions. The ex returned
610 * will already be evaluated. */
611 ex basic::subs(const exmap & m, unsigned options) const
616 // Substitute in subexpressions
617 for (size_t i=0; i<num; i++) {
618 const ex & orig_op = op(i);
619 const ex & subsed_op = orig_op.subs(m, options);
620 if (!are_ex_trivially_equal(orig_op, subsed_op)) {
622 // Something changed, clone the object
623 basic *copy = duplicate();
624 copy->setflag(status_flags::dynallocated);
625 copy->clearflag(status_flags::hash_calculated | status_flags::expanded);
627 // Substitute the changed operand
628 copy->let_op(i++) = subsed_op;
630 // Substitute the other operands
632 copy->let_op(i) = op(i).subs(m, options);
634 // Perform substitutions on the new object as a whole
635 return copy->subs_one_level(m, options);
640 // Nothing changed or no subexpressions
641 return subs_one_level(m, options);
644 /** Default interface of nth derivative ex::diff(s, n). It should be called
645 * instead of ::derivative(s) for first derivatives and for nth derivatives it
646 * just recurses down.
648 * @param s symbol to differentiate in
649 * @param nth order of differentiation
651 ex basic::diff(const symbol & s, unsigned nth) const
653 // trivial: zeroth derivative
657 // evaluate unevaluated *this before differentiating
658 if (!(flags & status_flags::evaluated))
659 return ex(*this).diff(s, nth);
661 ex ndiff = this->derivative(s);
662 while (!ndiff.is_zero() && // stop differentiating zeros
664 ndiff = ndiff.diff(s);
670 /** Return a vector containing the free indices of an expression. */
671 exvector basic::get_free_indices() const
673 return exvector(); // return an empty exvector
676 ex basic::eval_ncmul(const exvector & v) const
678 return hold_ncmul(v);
683 /** Function object to be applied by basic::derivative(). */
684 struct derivative_map_function : public map_function {
686 derivative_map_function(const symbol &sym) : s(sym) {}
687 ex operator()(const ex & e) { return diff(e, s); }
690 /** Default implementation of ex::diff(). It maps the operation on the
691 * operands (or returns 0 when the object has no operands).
694 ex basic::derivative(const symbol & s) const
699 derivative_map_function map_derivative(s);
700 return map(map_derivative);
704 /** Returns order relation between two objects of same type. This needs to be
705 * implemented by each class. It may never return anything else than 0,
706 * signalling equality, or +1 and -1 signalling inequality and determining
707 * the canonical ordering. (Perl hackers will wonder why C++ doesn't feature
708 * the spaceship operator <=> for denoting just this.) */
709 int basic::compare_same_type(const basic & other) const
711 return compare_pointers(this, &other);
714 /** Returns true if two objects of same type are equal. Normally needs
715 * not be reimplemented as long as it wasn't overwritten by some parent
716 * class, since it just calls compare_same_type(). The reason why this
717 * function exists is that sometimes it is easier to determine equality
718 * than an order relation and then it can be overridden. */
719 bool basic::is_equal_same_type(const basic & other) const
721 return compare_same_type(other)==0;
724 /** Returns true if the attributes of two objects are similar enough for
725 * a match. This function must not match subexpressions (this is already
726 * done by basic::match()). Only attributes not accessible by op() should
727 * be compared. This is also the reason why this function doesn't take the
728 * wildcard replacement list from match() as an argument: only subexpressions
729 * are subject to wildcard matches. Also, this function only needs to be
730 * implemented for container classes because is_equal_same_type() is
731 * automatically used instead of match_same_type() if nops() == 0.
733 * @see basic::match */
734 bool basic::match_same_type(const basic & other) const
736 // The default is to only consider subexpressions, but not any other
741 unsigned basic::return_type() const
743 return return_types::commutative;
746 unsigned basic::return_type_tinfo() const
751 /** Compute the hash value of an object and if it makes sense to store it in
752 * the objects status_flags, do so. The method inherited from class basic
753 * computes a hash value based on the type and hash values of possible
754 * members. For this reason it is well suited for container classes but
755 * atomic classes should override this implementation because otherwise they
756 * would all end up with the same hashvalue. */
757 unsigned basic::calchash() const
759 unsigned v = golden_ratio_hash(tinfo());
760 for (size_t i=0; i<nops(); i++) {
762 v ^= this->op(i).gethash();
765 // store calculated hash value only if object is already evaluated
766 if (flags & status_flags::evaluated) {
767 setflag(status_flags::hash_calculated);
774 /** Function object to be applied by basic::expand(). */
775 struct expand_map_function : public map_function {
777 expand_map_function(unsigned o) : options(o) {}
778 ex operator()(const ex & e) { return expand(e, options); }
781 /** Expand expression, i.e. multiply it out and return the result as a new
783 ex basic::expand(unsigned options) const
786 return (options == 0) ? setflag(status_flags::expanded) : *this;
788 expand_map_function map_expand(options);
789 return ex_to<basic>(map(map_expand)).setflag(options == 0 ? status_flags::expanded : 0);
795 // non-virtual functions in this class
800 /** Compare objects syntactically to establish canonical ordering.
801 * All compare functions return: -1 for *this less than other, 0 equal,
803 int basic::compare(const basic & other) const
805 const unsigned hash_this = gethash();
806 const unsigned hash_other = other.gethash();
807 if (hash_this<hash_other) return -1;
808 if (hash_this>hash_other) return 1;
810 const unsigned typeid_this = tinfo();
811 const unsigned typeid_other = other.tinfo();
812 if (typeid_this==typeid_other) {
813 GINAC_ASSERT(typeid(*this)==typeid(other));
814 // int cmpval = compare_same_type(other);
816 // std::cout << "hash collision, same type: "
817 // << *this << " and " << other << std::endl;
818 // this->print(print_tree(std::cout));
819 // std::cout << " and ";
820 // other.print(print_tree(std::cout));
821 // std::cout << std::endl;
824 return compare_same_type(other);
826 // std::cout << "hash collision, different types: "
827 // << *this << " and " << other << std::endl;
828 // this->print(print_tree(std::cout));
829 // std::cout << " and ";
830 // other.print(print_tree(std::cout));
831 // std::cout << std::endl;
832 return (typeid_this<typeid_other ? -1 : 1);
836 /** Test for syntactic equality.
837 * This is only a quick test, meaning objects should be in the same domain.
838 * You might have to .expand(), .normal() objects first, depending on the
839 * domain of your computation, to get a more reliable answer.
841 * @see is_equal_same_type */
842 bool basic::is_equal(const basic & other) const
844 if (this->gethash()!=other.gethash())
846 if (this->tinfo()!=other.tinfo())
849 GINAC_ASSERT(typeid(*this)==typeid(other));
851 return is_equal_same_type(other);
856 /** Stop further evaluation.
858 * @see basic::eval */
859 const basic & basic::hold() const
861 return setflag(status_flags::evaluated);
864 /** Ensure the object may be modified without hurting others, throws if this
865 * is not the case. */
866 void basic::ensure_if_modifiable() const
869 throw(std::runtime_error("cannot modify multiply referenced object"));
870 clearflag(status_flags::hash_calculated | status_flags::evaluated);
877 int max_recursion_level = 1024;