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 // Double dispatch on object type and print_context type
129 const registered_class_info * reg_info = &get_class_info();
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()
178 << std::hex << ", hash=0x" << hashvalue << ", flags=0x" << flags << std::dec
179 << ", 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(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) const
282 if (match(pattern, repl_lst))
284 for (size_t i=0; i<nops(); i++)
285 if (op(i).has(pattern))
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
299 basic *copy = duplicate();
300 copy->setflag(status_flags::dynallocated);
301 copy->clearflag(status_flags::hash_calculated | status_flags::expanded);
302 for (size_t i=0; i<num; i++)
303 copy->let_op(i) = f(copy->op(i));
307 /** Return degree of highest power in object s. */
308 int basic::degree(const ex & s) const
310 return is_equal(ex_to<basic>(s)) ? 1 : 0;
313 /** Return degree of lowest power in object s. */
314 int basic::ldegree(const ex & s) const
316 return is_equal(ex_to<basic>(s)) ? 1 : 0;
319 /** Return coefficient of degree n in object s. */
320 ex basic::coeff(const ex & s, int n) const
322 if (is_equal(ex_to<basic>(s)))
323 return n==1 ? _ex1 : _ex0;
325 return n==0 ? *this : _ex0;
328 /** Sort expanded expression in terms of powers of some object(s).
329 * @param s object(s) to sort in
330 * @param distributed recursive or distributed form (only used when s is a list) */
331 ex basic::collect(const ex & s, bool distributed) const
336 // List of objects specified
340 return collect(s.op(0));
342 else if (distributed) {
344 // Get lower/upper degree of all symbols in list
345 size_t num = s.nops();
349 int cnt; // current degree, 'counter'
350 ex coeff; // coefficient for degree 'cnt'
352 sym_info *si = new sym_info[num];
354 for (size_t i=0; i<num; i++) {
356 si[i].ldeg = si[i].cnt = this->ldegree(si[i].sym);
357 si[i].deg = this->degree(si[i].sym);
358 c = si[i].coeff = c.coeff(si[i].sym, si[i].cnt);
363 // Calculate coeff*x1^c1*...*xn^cn
365 for (size_t i=0; i<num; i++) {
367 y *= power(si[i].sym, cnt);
369 x += y * si[num - 1].coeff;
371 // Increment counters
375 if (si[n].cnt <= si[n].deg) {
376 // Update coefficients
382 for (size_t i=n; i<num; i++)
383 c = si[i].coeff = c.coeff(si[i].sym, si[i].cnt);
388 si[n].cnt = si[n].ldeg;
399 size_t n = s.nops() - 1;
410 // Only one object specified
411 for (int n=this->ldegree(s); n<=this->degree(s); ++n)
412 x += this->coeff(s,n)*power(s,n);
415 // correct for lost fractional arguments and return
416 return x + (*this - x).expand();
419 /** Perform automatic non-interruptive term rewriting rules. */
420 ex basic::eval(int level) const
422 // There is nothing to do for basic objects:
426 /** Function object to be applied by basic::evalf(). */
427 struct evalf_map_function : public map_function {
429 evalf_map_function(int l) : level(l) {}
430 ex operator()(const ex & e) { return evalf(e, level); }
433 /** Evaluate object numerically. */
434 ex basic::evalf(int level) const
441 else if (level == -max_recursion_level)
442 throw(std::runtime_error("max recursion level reached"));
444 evalf_map_function map_evalf(level - 1);
445 return map(map_evalf);
450 /** Function object to be applied by basic::evalm(). */
451 struct evalm_map_function : public map_function {
452 ex operator()(const ex & e) { return evalm(e); }
455 /** Evaluate sums, products and integer powers of matrices. */
456 ex basic::evalm() const
461 return map(map_evalm);
464 /** Perform automatic symbolic evaluations on indexed expression that
465 * contains this object as the base expression. */
466 ex basic::eval_indexed(const basic & i) const
467 // this function can't take a "const ex & i" because that would result
468 // in an infinite eval() loop
470 // There is nothing to do for basic objects
474 /** Add two indexed expressions. They are guaranteed to be of class indexed
475 * (or a subclass) and their indices are compatible. This function is used
476 * internally by simplify_indexed().
478 * @param self First indexed expression; it's base object is *this
479 * @param other Second indexed expression
480 * @return sum of self and other
481 * @see ex::simplify_indexed() */
482 ex basic::add_indexed(const ex & self, const ex & other) const
487 /** Multiply an indexed expression with a scalar. This function is used
488 * internally by simplify_indexed().
490 * @param self Indexed expression; it's base object is *this
491 * @param other Numeric value
492 * @return product of self and other
493 * @see ex::simplify_indexed() */
494 ex basic::scalar_mul_indexed(const ex & self, const numeric & other) const
499 /** Try to contract two indexed expressions that appear in the same product.
500 * If a contraction exists, the function overwrites one or both of the
501 * expressions and returns true. Otherwise it returns false. It is
502 * guaranteed that both expressions are of class indexed (or a subclass)
503 * and that at least one dummy index has been found. This functions is
504 * used internally by simplify_indexed().
506 * @param self Pointer to first indexed expression; it's base object is *this
507 * @param other Pointer to second indexed expression
508 * @param v The complete vector of factors
509 * @return true if the contraction was successful, false otherwise
510 * @see ex::simplify_indexed() */
511 bool basic::contract_with(exvector::iterator self, exvector::iterator other, exvector & v) const
517 /** Check whether the expression matches a given pattern. For every wildcard
518 * object in the pattern, an expression of the form "wildcard == matching_expression"
519 * is added to repl_lst. */
520 bool basic::match(const ex & pattern, lst & repl_lst) const
523 Sweet sweet shapes, sweet sweet shapes,
524 That's the key thing, right right.
525 Feed feed face, feed feed shapes,
526 But who is the king tonight?
527 Who is the king tonight?
528 Pattern is the thing, the key thing-a-ling,
529 But who is the king of Pattern?
530 But who is the king, the king thing-a-ling,
531 Who is the king of Pattern?
532 Bog is the king, the king thing-a-ling,
533 Bog is the king of Pattern.
534 Ba bu-bu-bu-bu bu-bu-bu-bu-bu-bu bu-bu
535 Bog is the king of Pattern.
538 if (is_exactly_a<wildcard>(pattern)) {
540 // Wildcard matches anything, but check whether we already have found
541 // a match for that wildcard first (if so, it the earlier match must
542 // be the same expression)
543 for (lst::const_iterator it = repl_lst.begin(); it != repl_lst.end(); ++it) {
544 if (it->op(0).is_equal(pattern))
545 return is_equal(ex_to<basic>(it->op(1)));
547 repl_lst.append(pattern == *this);
552 // Expression must be of the same type as the pattern
553 if (tinfo() != ex_to<basic>(pattern).tinfo())
556 // Number of subexpressions must match
557 if (nops() != pattern.nops())
560 // No subexpressions? Then just compare the objects (there can't be
561 // wildcards in the pattern)
563 return is_equal_same_type(ex_to<basic>(pattern));
565 // Check whether attributes that are not subexpressions match
566 if (!match_same_type(ex_to<basic>(pattern)))
569 // Otherwise the subexpressions must match one-to-one
570 for (size_t i=0; i<nops(); i++)
571 if (!op(i).match(pattern.op(i), repl_lst))
574 // Looks similar enough, match found
579 /** Helper function for subs(). Does not recurse into subexpressions. */
580 ex basic::subs_one_level(const exmap & m, unsigned options) const
582 exmap::const_iterator it;
584 if (options & subs_options::no_pattern) {
589 for (it = m.begin(); it != m.end(); ++it) {
591 if (match(ex_to<basic>(it->first), repl_lst))
592 return it->second.subs(repl_lst, options | subs_options::no_pattern); // avoid infinite recursion when re-substituting the wildcards
599 /** Substitute a set of objects by arbitrary expressions. The ex returned
600 * will already be evaluated. */
601 ex basic::subs(const exmap & m, unsigned options) const
606 // Substitute in subexpressions
607 for (size_t i=0; i<num; i++) {
608 const ex & orig_op = op(i);
609 const ex & subsed_op = orig_op.subs(m, options);
610 if (!are_ex_trivially_equal(orig_op, subsed_op)) {
612 // Something changed, clone the object
613 basic *copy = duplicate();
614 copy->setflag(status_flags::dynallocated);
615 copy->clearflag(status_flags::hash_calculated | status_flags::expanded);
617 // Substitute the changed operand
618 copy->let_op(i++) = subsed_op;
620 // Substitute the other operands
622 copy->let_op(i) = op(i).subs(m, options);
624 // Perform substitutions on the new object as a whole
625 return copy->subs_one_level(m, options);
630 // Nothing changed or no subexpressions
631 return subs_one_level(m, options);
634 /** Default interface of nth derivative ex::diff(s, n). It should be called
635 * instead of ::derivative(s) for first derivatives and for nth derivatives it
636 * just recurses down.
638 * @param s symbol to differentiate in
639 * @param nth order of differentiation
641 ex basic::diff(const symbol & s, unsigned nth) const
643 // trivial: zeroth derivative
647 // evaluate unevaluated *this before differentiating
648 if (!(flags & status_flags::evaluated))
649 return ex(*this).diff(s, nth);
651 ex ndiff = this->derivative(s);
652 while (!ndiff.is_zero() && // stop differentiating zeros
654 ndiff = ndiff.diff(s);
660 /** Return a vector containing the free indices of an expression. */
661 exvector basic::get_free_indices() const
663 return exvector(); // return an empty exvector
666 ex basic::eval_ncmul(const exvector & v) const
668 return hold_ncmul(v);
673 /** Function object to be applied by basic::derivative(). */
674 struct derivative_map_function : public map_function {
676 derivative_map_function(const symbol &sym) : s(sym) {}
677 ex operator()(const ex & e) { return diff(e, s); }
680 /** Default implementation of ex::diff(). It maps the operation on the
681 * operands (or returns 0 when the object has no operands).
684 ex basic::derivative(const symbol & s) const
689 derivative_map_function map_derivative(s);
690 return map(map_derivative);
694 /** Returns order relation between two objects of same type. This needs to be
695 * implemented by each class. It may never return anything else than 0,
696 * signalling equality, or +1 and -1 signalling inequality and determining
697 * the canonical ordering. (Perl hackers will wonder why C++ doesn't feature
698 * the spaceship operator <=> for denoting just this.) */
699 int basic::compare_same_type(const basic & other) const
701 return compare_pointers(this, &other);
704 /** Returns true if two objects of same type are equal. Normally needs
705 * not be reimplemented as long as it wasn't overwritten by some parent
706 * class, since it just calls compare_same_type(). The reason why this
707 * function exists is that sometimes it is easier to determine equality
708 * than an order relation and then it can be overridden. */
709 bool basic::is_equal_same_type(const basic & other) const
711 return compare_same_type(other)==0;
714 /** Returns true if the attributes of two objects are similar enough for
715 * a match. This function must not match subexpressions (this is already
716 * done by basic::match()). Only attributes not accessible by op() should
717 * be compared. This is also the reason why this function doesn't take the
718 * wildcard replacement list from match() as an argument: only subexpressions
719 * are subject to wildcard matches. Also, this function only needs to be
720 * implemented for container classes because is_equal_same_type() is
721 * automatically used instead of match_same_type() if nops() == 0.
723 * @see basic::match */
724 bool basic::match_same_type(const basic & other) const
726 // The default is to only consider subexpressions, but not any other
731 unsigned basic::return_type() const
733 return return_types::commutative;
736 unsigned basic::return_type_tinfo() const
741 /** Compute the hash value of an object and if it makes sense to store it in
742 * the objects status_flags, do so. The method inherited from class basic
743 * computes a hash value based on the type and hash values of possible
744 * members. For this reason it is well suited for container classes but
745 * atomic classes should override this implementation because otherwise they
746 * would all end up with the same hashvalue. */
747 unsigned basic::calchash() const
749 unsigned v = golden_ratio_hash(tinfo());
750 for (size_t i=0; i<nops(); i++) {
752 v ^= this->op(i).gethash();
755 // store calculated hash value only if object is already evaluated
756 if (flags & status_flags::evaluated) {
757 setflag(status_flags::hash_calculated);
764 /** Function object to be applied by basic::expand(). */
765 struct expand_map_function : public map_function {
767 expand_map_function(unsigned o) : options(o) {}
768 ex operator()(const ex & e) { return expand(e, options); }
771 /** Expand expression, i.e. multiply it out and return the result as a new
773 ex basic::expand(unsigned options) const
776 return (options == 0) ? setflag(status_flags::expanded) : *this;
778 expand_map_function map_expand(options);
779 return ex_to<basic>(map(map_expand)).setflag(options == 0 ? status_flags::expanded : 0);
785 // non-virtual functions in this class
790 /** Compare objects syntactically to establish canonical ordering.
791 * All compare functions return: -1 for *this less than other, 0 equal,
793 int basic::compare(const basic & other) const
795 const unsigned hash_this = gethash();
796 const unsigned hash_other = other.gethash();
797 if (hash_this<hash_other) return -1;
798 if (hash_this>hash_other) return 1;
800 const unsigned typeid_this = tinfo();
801 const unsigned typeid_other = other.tinfo();
802 if (typeid_this==typeid_other) {
803 GINAC_ASSERT(typeid(*this)==typeid(other));
804 // int cmpval = compare_same_type(other);
806 // std::cout << "hash collision, same type: "
807 // << *this << " and " << other << std::endl;
808 // this->print(print_tree(std::cout));
809 // std::cout << " and ";
810 // other.print(print_tree(std::cout));
811 // std::cout << std::endl;
814 return compare_same_type(other);
816 // std::cout << "hash collision, different types: "
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;
822 return (typeid_this<typeid_other ? -1 : 1);
826 /** Test for syntactic equality.
827 * This is only a quick test, meaning objects should be in the same domain.
828 * You might have to .expand(), .normal() objects first, depending on the
829 * domain of your computation, to get a more reliable answer.
831 * @see is_equal_same_type */
832 bool basic::is_equal(const basic & other) const
834 if (this->gethash()!=other.gethash())
836 if (this->tinfo()!=other.tinfo())
839 GINAC_ASSERT(typeid(*this)==typeid(other));
841 return is_equal_same_type(other);
846 /** Stop further evaluation.
848 * @see basic::eval */
849 const basic & basic::hold() const
851 return setflag(status_flags::evaluated);
854 /** Ensure the object may be modified without hurting others, throws if this
855 * is not the case. */
856 void basic::ensure_if_modifiable() const
859 throw(std::runtime_error("cannot modify multiply referenced object"));
860 clearflag(status_flags::hash_calculated | status_flags::evaluated);
867 int max_recursion_level = 1024;