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"
45 GINAC_IMPLEMENT_REGISTERED_CLASS(basic, void)
48 // default constructor, destructor, copy constructor and assignment operator
53 /** basic copy constructor: implicitly assumes that the other class is of
54 * the exact same type (as it's used by duplicate()), so it can copy the
55 * tinfo_key and the hash value. */
56 basic::basic(const basic & other) : tinfo_key(other.tinfo_key), flags(other.flags & ~status_flags::dynallocated), hashvalue(other.hashvalue), refcount(0)
58 GINAC_ASSERT(typeid(*this) == typeid(other));
61 /** basic assignment operator: the other object might be of a derived class. */
62 const basic & basic::operator=(const basic & other)
64 unsigned fl = other.flags & ~status_flags::dynallocated;
65 if (tinfo_key != other.tinfo_key) {
66 // The other object is of a derived class, so clear the flags as they
67 // might no longer apply (especially hash_calculated). Oh, and don't
68 // copy the tinfo_key: it is already set correctly for this object.
71 // The objects are of the exact same class, so copy the hash value.
72 hashvalue = other.hashvalue;
93 /** Construct object from archive_node. */
94 basic::basic(const archive_node &n, lst &sym_lst) : flags(0), refcount(0)
96 // Reconstruct tinfo_key from class name
97 std::string class_name;
98 if (n.find_string("class", class_name))
99 tinfo_key = find_tinfo_key(class_name);
101 throw (std::runtime_error("archive node contains no class name"));
104 /** Unarchive the object. */
105 DEFAULT_UNARCHIVE(basic)
107 /** Archive the object. */
108 void basic::archive(archive_node &n) const
110 n.add_string("class", class_name());
114 // new virtual functions which can be overridden by derived classes
119 /** Output to stream.
120 * @param c print context object that describes the output formatting
121 * @param level value that is used to identify the precedence or indentation
122 * level for placing parentheses and formatting */
123 void basic::print(const print_context & c, unsigned level) const
125 if (is_a<print_tree>(c)) {
127 c.s << std::string(level, ' ') << class_name()
128 << std::hex << ", hash=0x" << hashvalue << ", flags=0x" << flags << std::dec
129 << ", nops=" << nops()
131 for (size_t i=0; i<nops(); ++i)
132 op(i).print(c, level + static_cast<const print_tree &>(c).delta_indent);
135 c.s << "[" << class_name() << " object]";
138 /** Little wrapper around print to be called within a debugger.
139 * This is needed because you cannot call foo.print(cout) from within the
140 * debugger because it might not know what cout is. This method can be
141 * invoked with no argument and it will simply print to stdout.
144 * @see basic::dbgprinttree */
145 void basic::dbgprint() const
147 this->print(std::cerr);
148 std::cerr << std::endl;
151 /** Little wrapper around printtree to be called within a debugger.
153 * @see basic::dbgprint */
154 void basic::dbgprinttree() const
156 this->print(print_tree(std::cerr));
159 /** Return relative operator precedence (for parenthizing output). */
160 unsigned basic::precedence() const
165 /** Information about the object.
167 * @see class info_flags */
168 bool basic::info(unsigned inf) const
170 // all possible properties are false for basic objects
174 /** Number of operands/members. */
175 size_t basic::nops() const
177 // iterating from 0 to nops() on atomic objects should be an empty loop,
178 // and accessing their elements is a range error. Container objects should
183 /** Return operand/member at position i. */
184 ex basic::op(size_t i) const
186 throw(std::range_error(std::string("basic::op(): ") + class_name() + std::string(" has no operands")));
189 /** Return modifyable operand/member at position i. */
190 ex & basic::let_op(size_t i)
192 ensure_if_modifiable();
193 throw(std::range_error(std::string("basic::let_op(): ") + class_name() + std::string(" has no operands")));
196 ex basic::operator[](const ex & index) const
198 if (is_exactly_a<numeric>(index))
199 return op(static_cast<size_t>(ex_to<numeric>(index).to_int()));
201 throw(std::invalid_argument(std::string("non-numeric indices not supported by ") + class_name()));
204 ex basic::operator[](size_t i) const
209 ex & basic::operator[](const ex & index)
211 if (is_exactly_a<numeric>(index))
212 return let_op(ex_to<numeric>(index).to_int());
214 throw(std::invalid_argument(std::string("non-numeric indices not supported by ") + class_name()));
217 ex & basic::operator[](size_t i)
222 /** Test for occurrence of a pattern. An object 'has' a pattern if it matches
223 * the pattern itself or one of the children 'has' it. As a consequence
224 * (according to the definition of children) given e=x+y+z, e.has(x) is true
225 * but e.has(x+y) is false. */
226 bool basic::has(const ex & pattern) const
229 if (match(pattern, repl_lst))
231 for (size_t i=0; i<nops(); i++)
232 if (op(i).has(pattern))
238 /** Construct new expression by applying the specified function to all
239 * sub-expressions (one level only, not recursively). */
240 ex basic::map(map_function & f) const
246 basic *copy = duplicate();
247 copy->setflag(status_flags::dynallocated);
248 copy->clearflag(status_flags::hash_calculated | status_flags::expanded);
249 for (size_t i=0; i<num; i++)
250 copy->let_op(i) = f(copy->op(i));
254 /** Return degree of highest power in object s. */
255 int basic::degree(const ex & s) const
257 return is_equal(ex_to<basic>(s)) ? 1 : 0;
260 /** Return degree of lowest power in object s. */
261 int basic::ldegree(const ex & s) const
263 return is_equal(ex_to<basic>(s)) ? 1 : 0;
266 /** Return coefficient of degree n in object s. */
267 ex basic::coeff(const ex & s, int n) const
269 if (is_equal(ex_to<basic>(s)))
270 return n==1 ? _ex1 : _ex0;
272 return n==0 ? *this : _ex0;
275 /** Sort expanded expression in terms of powers of some object(s).
276 * @param s object(s) to sort in
277 * @param distributed recursive or distributed form (only used when s is a list) */
278 ex basic::collect(const ex & s, bool distributed) const
283 // List of objects specified
287 return collect(s.op(0));
289 else if (distributed) {
291 // Get lower/upper degree of all symbols in list
292 size_t num = s.nops();
296 int cnt; // current degree, 'counter'
297 ex coeff; // coefficient for degree 'cnt'
299 sym_info *si = new sym_info[num];
301 for (size_t i=0; i<num; i++) {
303 si[i].ldeg = si[i].cnt = this->ldegree(si[i].sym);
304 si[i].deg = this->degree(si[i].sym);
305 c = si[i].coeff = c.coeff(si[i].sym, si[i].cnt);
310 // Calculate coeff*x1^c1*...*xn^cn
312 for (size_t i=0; i<num; i++) {
314 y *= power(si[i].sym, cnt);
316 x += y * si[num - 1].coeff;
318 // Increment counters
322 if (si[n].cnt <= si[n].deg) {
323 // Update coefficients
329 for (size_t i=n; i<num; i++)
330 c = si[i].coeff = c.coeff(si[i].sym, si[i].cnt);
335 si[n].cnt = si[n].ldeg;
346 size_t n = s.nops() - 1;
357 // Only one object specified
358 for (int n=this->ldegree(s); n<=this->degree(s); ++n)
359 x += this->coeff(s,n)*power(s,n);
362 // correct for lost fractional arguments and return
363 return x + (*this - x).expand();
366 /** Perform automatic non-interruptive term rewriting rules. */
367 ex basic::eval(int level) const
369 // There is nothing to do for basic objects:
373 /** Function object to be applied by basic::evalf(). */
374 struct evalf_map_function : public map_function {
376 evalf_map_function(int l) : level(l) {}
377 ex operator()(const ex & e) { return evalf(e, level); }
380 /** Evaluate object numerically. */
381 ex basic::evalf(int level) const
388 else if (level == -max_recursion_level)
389 throw(std::runtime_error("max recursion level reached"));
391 evalf_map_function map_evalf(level - 1);
392 return map(map_evalf);
397 /** Function object to be applied by basic::evalm(). */
398 struct evalm_map_function : public map_function {
399 ex operator()(const ex & e) { return evalm(e); }
402 /** Evaluate sums, products and integer powers of matrices. */
403 ex basic::evalm() const
408 return map(map_evalm);
411 /** Perform automatic symbolic evaluations on indexed expression that
412 * contains this object as the base expression. */
413 ex basic::eval_indexed(const basic & i) const
414 // this function can't take a "const ex & i" because that would result
415 // in an infinite eval() loop
417 // There is nothing to do for basic objects
421 /** Add two indexed expressions. They are guaranteed to be of class indexed
422 * (or a subclass) and their indices are compatible. This function is used
423 * internally by simplify_indexed().
425 * @param self First indexed expression; it's base object is *this
426 * @param other Second indexed expression
427 * @return sum of self and other
428 * @see ex::simplify_indexed() */
429 ex basic::add_indexed(const ex & self, const ex & other) const
434 /** Multiply an indexed expression with a scalar. This function is used
435 * internally by simplify_indexed().
437 * @param self Indexed expression; it's base object is *this
438 * @param other Numeric value
439 * @return product of self and other
440 * @see ex::simplify_indexed() */
441 ex basic::scalar_mul_indexed(const ex & self, const numeric & other) const
446 /** Try to contract two indexed expressions that appear in the same product.
447 * If a contraction exists, the function overwrites one or both of the
448 * expressions and returns true. Otherwise it returns false. It is
449 * guaranteed that both expressions are of class indexed (or a subclass)
450 * and that at least one dummy index has been found. This functions is
451 * used internally by simplify_indexed().
453 * @param self Pointer to first indexed expression; it's base object is *this
454 * @param other Pointer to second indexed expression
455 * @param v The complete vector of factors
456 * @return true if the contraction was successful, false otherwise
457 * @see ex::simplify_indexed() */
458 bool basic::contract_with(exvector::iterator self, exvector::iterator other, exvector & v) const
464 /** Check whether the expression matches a given pattern. For every wildcard
465 * object in the pattern, an expression of the form "wildcard == matching_expression"
466 * is added to repl_lst. */
467 bool basic::match(const ex & pattern, lst & repl_lst) const
470 Sweet sweet shapes, sweet sweet shapes,
471 That's the key thing, right right.
472 Feed feed face, feed feed shapes,
473 But who is the king tonight?
474 Who is the king tonight?
475 Pattern is the thing, the key thing-a-ling,
476 But who is the king of Pattern?
477 But who is the king, the king thing-a-ling,
478 Who is the king of Pattern?
479 Bog is the king, the king thing-a-ling,
480 Bog is the king of Pattern.
481 Ba bu-bu-bu-bu bu-bu-bu-bu-bu-bu bu-bu
482 Bog is the king of Pattern.
485 if (is_exactly_a<wildcard>(pattern)) {
487 // Wildcard matches anything, but check whether we already have found
488 // a match for that wildcard first (if so, it the earlier match must
489 // be the same expression)
490 for (lst::const_iterator it = repl_lst.begin(); it != repl_lst.end(); ++it) {
491 if (it->op(0).is_equal(pattern))
492 return is_equal(ex_to<basic>(it->op(1)));
494 repl_lst.append(pattern == *this);
499 // Expression must be of the same type as the pattern
500 if (tinfo() != ex_to<basic>(pattern).tinfo())
503 // Number of subexpressions must match
504 if (nops() != pattern.nops())
507 // No subexpressions? Then just compare the objects (there can't be
508 // wildcards in the pattern)
510 return is_equal_same_type(ex_to<basic>(pattern));
512 // Check whether attributes that are not subexpressions match
513 if (!match_same_type(ex_to<basic>(pattern)))
516 // Otherwise the subexpressions must match one-to-one
517 for (size_t i=0; i<nops(); i++)
518 if (!op(i).match(pattern.op(i), repl_lst))
521 // Looks similar enough, match found
526 /** Helper function for subs(). Does not recurse into subexpressions. */
527 ex basic::subs_one_level(const lst & ls, const lst & lr, unsigned options) const
529 GINAC_ASSERT(ls.nops() == lr.nops());
531 lst::const_iterator its, itr;
533 if (options & subs_options::subs_no_pattern) {
534 for (its = ls.begin(), itr = lr.begin(); its != ls.end(); ++its, ++itr) {
535 if (is_equal(ex_to<basic>(*its)))
539 for (its = ls.begin(), itr = lr.begin(); its != ls.end(); ++its, ++itr) {
541 if (match(ex_to<basic>(*its), repl_lst))
542 return itr->subs(repl_lst, options | subs_options::subs_no_pattern); // avoid infinite recursion when re-substituting the wildcards
549 /** Substitute a set of objects by arbitrary expressions. The ex returned
550 * will already be evaluated. */
551 ex basic::subs(const lst & ls, const lst & lr, unsigned options) const
556 // Substitute in subexpressions
557 for (size_t i=0; i<num; i++) {
558 const ex & orig_op = op(i);
559 const ex & subsed_op = orig_op.subs(ls, lr, options);
560 if (!are_ex_trivially_equal(orig_op, subsed_op)) {
562 // Something changed, clone the object
563 basic *copy = duplicate();
564 copy->setflag(status_flags::dynallocated);
565 copy->clearflag(status_flags::hash_calculated | status_flags::expanded);
567 // Substitute the changed operand
568 copy->let_op(i++) = subsed_op;
570 // Substitute the other operands
572 copy->let_op(i) = op(i).subs(ls, lr, options);
574 // Perform substitutions on the new object as a whole
575 return copy->subs_one_level(ls, lr, options);
580 // Nothing changed or no subexpressions
581 return subs_one_level(ls, lr, options);
584 /** Default interface of nth derivative ex::diff(s, n). It should be called
585 * instead of ::derivative(s) for first derivatives and for nth derivatives it
586 * just recurses down.
588 * @param s symbol to differentiate in
589 * @param nth order of differentiation
591 ex basic::diff(const symbol & s, unsigned nth) const
593 // trivial: zeroth derivative
597 // evaluate unevaluated *this before differentiating
598 if (!(flags & status_flags::evaluated))
599 return ex(*this).diff(s, nth);
601 ex ndiff = this->derivative(s);
602 while (!ndiff.is_zero() && // stop differentiating zeros
604 ndiff = ndiff.diff(s);
610 /** Return a vector containing the free indices of an expression. */
611 exvector basic::get_free_indices() const
613 return exvector(); // return an empty exvector
616 ex basic::eval_ncmul(const exvector & v) const
618 return hold_ncmul(v);
623 /** Function object to be applied by basic::derivative(). */
624 struct derivative_map_function : public map_function {
626 derivative_map_function(const symbol &sym) : s(sym) {}
627 ex operator()(const ex & e) { return diff(e, s); }
630 /** Default implementation of ex::diff(). It maps the operation on the
631 * operands (or returns 0 when the object has no operands).
634 ex basic::derivative(const symbol & s) const
639 derivative_map_function map_derivative(s);
640 return map(map_derivative);
644 /** Returns order relation between two objects of same type. This needs to be
645 * implemented by each class. It may never return anything else than 0,
646 * signalling equality, or +1 and -1 signalling inequality and determining
647 * the canonical ordering. (Perl hackers will wonder why C++ doesn't feature
648 * the spaceship operator <=> for denoting just this.) */
649 int basic::compare_same_type(const basic & other) const
651 return compare_pointers(this, &other);
654 /** Returns true if two objects of same type are equal. Normally needs
655 * not be reimplemented as long as it wasn't overwritten by some parent
656 * class, since it just calls compare_same_type(). The reason why this
657 * function exists is that sometimes it is easier to determine equality
658 * than an order relation and then it can be overridden. */
659 bool basic::is_equal_same_type(const basic & other) const
661 return compare_same_type(other)==0;
664 /** Returns true if the attributes of two objects are similar enough for
665 * a match. This function must not match subexpressions (this is already
666 * done by basic::match()). Only attributes not accessible by op() should
667 * be compared. This is also the reason why this function doesn't take the
668 * wildcard replacement list from match() as an argument: only subexpressions
669 * are subject to wildcard matches. Also, this function only needs to be
670 * implemented for container classes because is_equal_same_type() is
671 * automatically used instead of match_same_type() if nops() == 0.
673 * @see basic::match */
674 bool basic::match_same_type(const basic & other) const
676 // The default is to only consider subexpressions, but not any other
681 unsigned basic::return_type() const
683 return return_types::commutative;
686 unsigned basic::return_type_tinfo() const
691 /** Compute the hash value of an object and if it makes sense to store it in
692 * the objects status_flags, do so. The method inherited from class basic
693 * computes a hash value based on the type and hash values of possible
694 * members. For this reason it is well suited for container classes but
695 * atomic classes should override this implementation because otherwise they
696 * would all end up with the same hashvalue. */
697 unsigned basic::calchash() const
699 unsigned v = golden_ratio_hash(tinfo());
700 for (size_t i=0; i<nops(); i++) {
702 v ^= (const_cast<basic *>(this))->op(i).gethash();
705 // store calculated hash value only if object is already evaluated
706 if (flags & status_flags::evaluated) {
707 setflag(status_flags::hash_calculated);
714 /** Function object to be applied by basic::expand(). */
715 struct expand_map_function : public map_function {
717 expand_map_function(unsigned o) : options(o) {}
718 ex operator()(const ex & e) { return expand(e, options); }
721 /** Expand expression, i.e. multiply it out and return the result as a new
723 ex basic::expand(unsigned options) const
726 return (options == 0) ? setflag(status_flags::expanded) : *this;
728 expand_map_function map_expand(options);
729 return ex_to<basic>(map(map_expand)).setflag(options == 0 ? status_flags::expanded : 0);
735 // non-virtual functions in this class
740 /** Substitute objects in an expression (syntactic substitution) and return
741 * the result as a new expression. There are two valid types of
742 * replacement arguments: 1) a relational like object==ex and 2) a list of
743 * relationals lst(object1==ex1,object2==ex2,...), which is converted to
744 * subs(lst(object1,object2,...),lst(ex1,ex2,...)). */
745 ex basic::subs(const ex & e, unsigned options) const
747 if (e.info(info_flags::relation_equal)) {
748 return subs(lst(e), options);
750 if (!e.info(info_flags::list)) {
751 throw(std::invalid_argument("basic::subs(ex): argument must be a list"));
754 // Split list into two
757 GINAC_ASSERT(is_a<lst>(e));
758 for (lst::const_iterator it = ex_to<lst>(e).begin(); it != ex_to<lst>(e).end(); ++it) {
760 if (!r.info(info_flags::relation_equal)) {
761 throw(std::invalid_argument("basic::subs(ex): argument must be a list of equations"));
766 return subs(ls, lr, options);
769 /** Compare objects syntactically to establish canonical ordering.
770 * All compare functions return: -1 for *this less than other, 0 equal,
772 int basic::compare(const basic & other) const
774 const unsigned hash_this = gethash();
775 const unsigned hash_other = other.gethash();
776 if (hash_this<hash_other) return -1;
777 if (hash_this>hash_other) return 1;
779 const unsigned typeid_this = tinfo();
780 const unsigned typeid_other = other.tinfo();
781 if (typeid_this==typeid_other) {
782 GINAC_ASSERT(typeid(*this)==typeid(other));
783 // int cmpval = compare_same_type(other);
785 // std::cout << "hash collision, same type: "
786 // << *this << " and " << other << std::endl;
787 // this->print(print_tree(std::cout));
788 // std::cout << " and ";
789 // other.print(print_tree(std::cout));
790 // std::cout << std::endl;
793 return compare_same_type(other);
795 // std::cout << "hash collision, different types: "
796 // << *this << " and " << other << std::endl;
797 // this->print(print_tree(std::cout));
798 // std::cout << " and ";
799 // other.print(print_tree(std::cout));
800 // std::cout << std::endl;
801 return (typeid_this<typeid_other ? -1 : 1);
805 /** Test for syntactic equality.
806 * This is only a quick test, meaning objects should be in the same domain.
807 * You might have to .expand(), .normal() objects first, depending on the
808 * domain of your computation, to get a more reliable answer.
810 * @see is_equal_same_type */
811 bool basic::is_equal(const basic & other) const
813 if (this->gethash()!=other.gethash())
815 if (this->tinfo()!=other.tinfo())
818 GINAC_ASSERT(typeid(*this)==typeid(other));
820 return is_equal_same_type(other);
825 /** Stop further evaluation.
827 * @see basic::eval */
828 const basic & basic::hold() const
830 return setflag(status_flags::evaluated);
833 /** Ensure the object may be modified without hurting others, throws if this
834 * is not the case. */
835 void basic::ensure_if_modifiable() const
838 throw(std::runtime_error("cannot modify multiply referenced object"));
839 clearflag(status_flags::hash_calculated | status_flags::evaluated);
846 int max_recursion_level = 1024;