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 std::clog << "searching class " << reg_info->options.get_name() << std::endl;
144 const std::vector<print_functor> & pdt = reg_info->options.get_print_dispatch_table();
147 std::clog << "searching context " << pc_info->options.get_name() << ", ID " << pc_info->options.get_id() << std::endl;
148 unsigned id = pc_info->options.get_id();
149 if (id >= pdt.size() || !(pdt[id].is_valid())) {
151 // Method not found, try parent print_context class
152 const print_context_class_info * parent_pc_info = pc_info->get_parent();
153 if (parent_pc_info) {
154 pc_info = parent_pc_info;
158 // Method still not found, try parent class
159 const registered_class_info * parent_reg_info = reg_info->get_parent();
160 if (parent_reg_info) {
161 reg_info = parent_reg_info;
162 pc_info = &c.get_class_info();
166 // Method still not found. This shouldn't happen because basic (the
167 // base class of the algebraic hierarchy) registers a method for
168 // print_context (the base class of the print context hierarchy),
169 // so if we end up here, there's something wrong with the class
171 throw (std::runtime_error(std::string("basic::print(): method for ") + class_name() + "/" + c.class_name() + " not found"));
176 std::clog << "method found, calling" << std::endl;
177 std::clog << " this = " << class_name() << ", context = " << c.class_name() << std::endl;
178 pdt[id](*this, c, level);
182 /** Default output to stream. */
183 void basic::do_print(const print_context & c, unsigned level) const
185 c.s << "[" << class_name() << " object]";
188 /** Tree output to stream. */
189 void basic::do_print_tree(const print_tree & c, unsigned level) const
191 c.s << std::string(level, ' ') << class_name()
192 << std::hex << ", hash=0x" << hashvalue << ", flags=0x" << flags << std::dec
193 << ", nops=" << nops()
195 for (size_t i=0; i<nops(); ++i)
196 op(i).print(c, level + c.delta_indent);
199 /** Python parsable output to stream. */
200 void basic::do_print_python_repr(const print_python_repr & c, unsigned level) const
202 c.s << class_name() << "()";
205 /** Little wrapper around print to be called within a debugger.
206 * This is needed because you cannot call foo.print(cout) from within the
207 * debugger because it might not know what cout is. This method can be
208 * invoked with no argument and it will simply print to stdout.
211 * @see basic::dbgprinttree */
212 void basic::dbgprint() const
214 this->print(std::cerr);
215 std::cerr << std::endl;
218 /** Little wrapper around printtree to be called within a debugger.
220 * @see basic::dbgprint */
221 void basic::dbgprinttree() const
223 this->print(print_tree(std::cerr));
226 /** Return relative operator precedence (for parenthezing output). */
227 unsigned basic::precedence() const
232 /** Information about the object.
234 * @see class info_flags */
235 bool basic::info(unsigned inf) const
237 // all possible properties are false for basic objects
241 /** Number of operands/members. */
242 size_t basic::nops() const
244 // iterating from 0 to nops() on atomic objects should be an empty loop,
245 // and accessing their elements is a range error. Container objects should
250 /** Return operand/member at position i. */
251 ex basic::op(size_t i) const
253 throw(std::range_error(std::string("basic::op(): ") + class_name() + std::string(" has no operands")));
256 /** Return modifyable operand/member at position i. */
257 ex & basic::let_op(size_t i)
259 ensure_if_modifiable();
260 throw(std::range_error(std::string("basic::let_op(): ") + class_name() + std::string(" has no operands")));
263 ex basic::operator[](const ex & index) const
265 if (is_exactly_a<numeric>(index))
266 return op(static_cast<size_t>(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) const
276 ex & basic::operator[](const ex & index)
278 if (is_exactly_a<numeric>(index))
279 return let_op(ex_to<numeric>(index).to_int());
281 throw(std::invalid_argument(std::string("non-numeric indices not supported by ") + class_name()));
284 ex & basic::operator[](size_t i)
289 /** Test for occurrence of a pattern. An object 'has' a pattern if it matches
290 * the pattern itself or one of the children 'has' it. As a consequence
291 * (according to the definition of children) given e=x+y+z, e.has(x) is true
292 * but e.has(x+y) is false. */
293 bool basic::has(const ex & pattern) const
296 if (match(pattern, repl_lst))
298 for (size_t i=0; i<nops(); i++)
299 if (op(i).has(pattern))
305 /** Construct new expression by applying the specified function to all
306 * sub-expressions (one level only, not recursively). */
307 ex basic::map(map_function & f) const
313 basic *copy = duplicate();
314 copy->setflag(status_flags::dynallocated);
315 copy->clearflag(status_flags::hash_calculated | status_flags::expanded);
316 for (size_t i=0; i<num; i++)
317 copy->let_op(i) = f(copy->op(i));
321 /** Return degree of highest power in object s. */
322 int basic::degree(const ex & s) const
324 return is_equal(ex_to<basic>(s)) ? 1 : 0;
327 /** Return degree of lowest power in object s. */
328 int basic::ldegree(const ex & s) const
330 return is_equal(ex_to<basic>(s)) ? 1 : 0;
333 /** Return coefficient of degree n in object s. */
334 ex basic::coeff(const ex & s, int n) const
336 if (is_equal(ex_to<basic>(s)))
337 return n==1 ? _ex1 : _ex0;
339 return n==0 ? *this : _ex0;
342 /** Sort expanded expression in terms of powers of some object(s).
343 * @param s object(s) to sort in
344 * @param distributed recursive or distributed form (only used when s is a list) */
345 ex basic::collect(const ex & s, bool distributed) const
350 // List of objects specified
354 return collect(s.op(0));
356 else if (distributed) {
358 // Get lower/upper degree of all symbols in list
359 size_t num = s.nops();
363 int cnt; // current degree, 'counter'
364 ex coeff; // coefficient for degree 'cnt'
366 sym_info *si = new sym_info[num];
368 for (size_t i=0; i<num; i++) {
370 si[i].ldeg = si[i].cnt = this->ldegree(si[i].sym);
371 si[i].deg = this->degree(si[i].sym);
372 c = si[i].coeff = c.coeff(si[i].sym, si[i].cnt);
377 // Calculate coeff*x1^c1*...*xn^cn
379 for (size_t i=0; i<num; i++) {
381 y *= power(si[i].sym, cnt);
383 x += y * si[num - 1].coeff;
385 // Increment counters
389 if (si[n].cnt <= si[n].deg) {
390 // Update coefficients
396 for (size_t i=n; i<num; i++)
397 c = si[i].coeff = c.coeff(si[i].sym, si[i].cnt);
402 si[n].cnt = si[n].ldeg;
413 size_t n = s.nops() - 1;
424 // Only one object specified
425 for (int n=this->ldegree(s); n<=this->degree(s); ++n)
426 x += this->coeff(s,n)*power(s,n);
429 // correct for lost fractional arguments and return
430 return x + (*this - x).expand();
433 /** Perform automatic non-interruptive term rewriting rules. */
434 ex basic::eval(int level) const
436 // There is nothing to do for basic objects:
440 /** Function object to be applied by basic::evalf(). */
441 struct evalf_map_function : public map_function {
443 evalf_map_function(int l) : level(l) {}
444 ex operator()(const ex & e) { return evalf(e, level); }
447 /** Evaluate object numerically. */
448 ex basic::evalf(int level) const
455 else if (level == -max_recursion_level)
456 throw(std::runtime_error("max recursion level reached"));
458 evalf_map_function map_evalf(level - 1);
459 return map(map_evalf);
464 /** Function object to be applied by basic::evalm(). */
465 struct evalm_map_function : public map_function {
466 ex operator()(const ex & e) { return evalm(e); }
469 /** Evaluate sums, products and integer powers of matrices. */
470 ex basic::evalm() const
475 return map(map_evalm);
478 /** Perform automatic symbolic evaluations on indexed expression that
479 * contains this object as the base expression. */
480 ex basic::eval_indexed(const basic & i) const
481 // this function can't take a "const ex & i" because that would result
482 // in an infinite eval() loop
484 // There is nothing to do for basic objects
488 /** Add two indexed expressions. They are guaranteed to be of class indexed
489 * (or a subclass) and their indices are compatible. This function is used
490 * internally by simplify_indexed().
492 * @param self First indexed expression; it's base object is *this
493 * @param other Second indexed expression
494 * @return sum of self and other
495 * @see ex::simplify_indexed() */
496 ex basic::add_indexed(const ex & self, const ex & other) const
501 /** Multiply an indexed expression with a scalar. This function is used
502 * internally by simplify_indexed().
504 * @param self Indexed expression; it's base object is *this
505 * @param other Numeric value
506 * @return product of self and other
507 * @see ex::simplify_indexed() */
508 ex basic::scalar_mul_indexed(const ex & self, const numeric & other) const
513 /** Try to contract two indexed expressions that appear in the same product.
514 * If a contraction exists, the function overwrites one or both of the
515 * expressions and returns true. Otherwise it returns false. It is
516 * guaranteed that both expressions are of class indexed (or a subclass)
517 * and that at least one dummy index has been found. This functions is
518 * used internally by simplify_indexed().
520 * @param self Pointer to first indexed expression; it's base object is *this
521 * @param other Pointer to second indexed expression
522 * @param v The complete vector of factors
523 * @return true if the contraction was successful, false otherwise
524 * @see ex::simplify_indexed() */
525 bool basic::contract_with(exvector::iterator self, exvector::iterator other, exvector & v) const
531 /** Check whether the expression matches a given pattern. For every wildcard
532 * object in the pattern, an expression of the form "wildcard == matching_expression"
533 * is added to repl_lst. */
534 bool basic::match(const ex & pattern, lst & repl_lst) const
537 Sweet sweet shapes, sweet sweet shapes,
538 That's the key thing, right right.
539 Feed feed face, feed feed shapes,
540 But who is the king tonight?
541 Who is the king tonight?
542 Pattern is the thing, the key thing-a-ling,
543 But who is the king of Pattern?
544 But who is the king, the king thing-a-ling,
545 Who is the king of Pattern?
546 Bog is the king, the king thing-a-ling,
547 Bog is the king of Pattern.
548 Ba bu-bu-bu-bu bu-bu-bu-bu-bu-bu bu-bu
549 Bog is the king of Pattern.
552 if (is_exactly_a<wildcard>(pattern)) {
554 // Wildcard matches anything, but check whether we already have found
555 // a match for that wildcard first (if so, it the earlier match must
556 // be the same expression)
557 for (lst::const_iterator it = repl_lst.begin(); it != repl_lst.end(); ++it) {
558 if (it->op(0).is_equal(pattern))
559 return is_equal(ex_to<basic>(it->op(1)));
561 repl_lst.append(pattern == *this);
566 // Expression must be of the same type as the pattern
567 if (tinfo() != ex_to<basic>(pattern).tinfo())
570 // Number of subexpressions must match
571 if (nops() != pattern.nops())
574 // No subexpressions? Then just compare the objects (there can't be
575 // wildcards in the pattern)
577 return is_equal_same_type(ex_to<basic>(pattern));
579 // Check whether attributes that are not subexpressions match
580 if (!match_same_type(ex_to<basic>(pattern)))
583 // Otherwise the subexpressions must match one-to-one
584 for (size_t i=0; i<nops(); i++)
585 if (!op(i).match(pattern.op(i), repl_lst))
588 // 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 exmap::const_iterator it;
598 if (options & subs_options::no_pattern) {
603 for (it = m.begin(); it != m.end(); ++it) {
605 if (match(ex_to<basic>(it->first), repl_lst))
606 return it->second.subs(repl_lst, options | subs_options::no_pattern); // 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::eval_ncmul(const exvector & v) const
682 return hold_ncmul(v);
687 /** Function object to be applied by basic::derivative(). */
688 struct derivative_map_function : public map_function {
690 derivative_map_function(const symbol &sym) : s(sym) {}
691 ex operator()(const ex & e) { return diff(e, s); }
694 /** Default implementation of ex::diff(). It maps the operation on the
695 * operands (or returns 0 when the object has no operands).
698 ex basic::derivative(const symbol & s) const
703 derivative_map_function map_derivative(s);
704 return map(map_derivative);
708 /** Returns order relation between two objects of same type. This needs to be
709 * implemented by each class. It may never return anything else than 0,
710 * signalling equality, or +1 and -1 signalling inequality and determining
711 * the canonical ordering. (Perl hackers will wonder why C++ doesn't feature
712 * the spaceship operator <=> for denoting just this.) */
713 int basic::compare_same_type(const basic & other) const
715 return compare_pointers(this, &other);
718 /** Returns true if two objects of same type are equal. Normally needs
719 * not be reimplemented as long as it wasn't overwritten by some parent
720 * class, since it just calls compare_same_type(). The reason why this
721 * function exists is that sometimes it is easier to determine equality
722 * than an order relation and then it can be overridden. */
723 bool basic::is_equal_same_type(const basic & other) const
725 return compare_same_type(other)==0;
728 /** Returns true if the attributes of two objects are similar enough for
729 * a match. This function must not match subexpressions (this is already
730 * done by basic::match()). Only attributes not accessible by op() should
731 * be compared. This is also the reason why this function doesn't take the
732 * wildcard replacement list from match() as an argument: only subexpressions
733 * are subject to wildcard matches. Also, this function only needs to be
734 * implemented for container classes because is_equal_same_type() is
735 * automatically used instead of match_same_type() if nops() == 0.
737 * @see basic::match */
738 bool basic::match_same_type(const basic & other) const
740 // The default is to only consider subexpressions, but not any other
745 unsigned basic::return_type() const
747 return return_types::commutative;
750 unsigned basic::return_type_tinfo() const
755 /** Compute the hash value of an object and if it makes sense to store it in
756 * the objects status_flags, do so. The method inherited from class basic
757 * computes a hash value based on the type and hash values of possible
758 * members. For this reason it is well suited for container classes but
759 * atomic classes should override this implementation because otherwise they
760 * would all end up with the same hashvalue. */
761 unsigned basic::calchash() const
763 unsigned v = golden_ratio_hash(tinfo());
764 for (size_t i=0; i<nops(); i++) {
766 v ^= this->op(i).gethash();
769 // store calculated hash value only if object is already evaluated
770 if (flags & status_flags::evaluated) {
771 setflag(status_flags::hash_calculated);
778 /** Function object to be applied by basic::expand(). */
779 struct expand_map_function : public map_function {
781 expand_map_function(unsigned o) : options(o) {}
782 ex operator()(const ex & e) { return expand(e, options); }
785 /** Expand expression, i.e. multiply it out and return the result as a new
787 ex basic::expand(unsigned options) const
790 return (options == 0) ? setflag(status_flags::expanded) : *this;
792 expand_map_function map_expand(options);
793 return ex_to<basic>(map(map_expand)).setflag(options == 0 ? status_flags::expanded : 0);
799 // non-virtual functions in this class
804 /** Compare objects syntactically to establish canonical ordering.
805 * All compare functions return: -1 for *this less than other, 0 equal,
807 int basic::compare(const basic & other) const
809 const unsigned hash_this = gethash();
810 const unsigned hash_other = other.gethash();
811 if (hash_this<hash_other) return -1;
812 if (hash_this>hash_other) return 1;
814 const unsigned typeid_this = tinfo();
815 const unsigned typeid_other = other.tinfo();
816 if (typeid_this==typeid_other) {
817 GINAC_ASSERT(typeid(*this)==typeid(other));
818 // int cmpval = compare_same_type(other);
820 // std::cout << "hash collision, same type: "
821 // << *this << " and " << other << std::endl;
822 // this->print(print_tree(std::cout));
823 // std::cout << " and ";
824 // other.print(print_tree(std::cout));
825 // std::cout << std::endl;
828 return compare_same_type(other);
830 // std::cout << "hash collision, different types: "
831 // << *this << " and " << other << std::endl;
832 // this->print(print_tree(std::cout));
833 // std::cout << " and ";
834 // other.print(print_tree(std::cout));
835 // std::cout << std::endl;
836 return (typeid_this<typeid_other ? -1 : 1);
840 /** Test for syntactic equality.
841 * This is only a quick test, meaning objects should be in the same domain.
842 * You might have to .expand(), .normal() objects first, depending on the
843 * domain of your computation, to get a more reliable answer.
845 * @see is_equal_same_type */
846 bool basic::is_equal(const basic & other) const
848 if (this->gethash()!=other.gethash())
850 if (this->tinfo()!=other.tinfo())
853 GINAC_ASSERT(typeid(*this)==typeid(other));
855 return is_equal_same_type(other);
860 /** Stop further evaluation.
862 * @see basic::eval */
863 const basic & basic::hold() const
865 return setflag(status_flags::evaluated);
868 /** Ensure the object may be modified without hurting others, throws if this
869 * is not the case. */
870 void basic::ensure_if_modifiable() const
873 throw(std::runtime_error("cannot modify multiply referenced object"));
874 clearflag(status_flags::hash_calculated | status_flags::evaluated);
881 int max_recursion_level = 1024;