3 * Implementation of GiNaC's ABC. */
6 * GiNaC Copyright (C) 1999-2002 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"
44 GINAC_IMPLEMENT_REGISTERED_CLASS_NO_CTORS(basic, void)
47 // default ctor, dtor, copy ctor, assignment operator and helpers
52 basic::basic(const basic & other) : tinfo_key(TINFO_basic), flags(0), refcount(0)
57 const basic & basic::operator=(const basic & other)
68 // none (all conditionally inlined)
74 // none (all conditionally inlined)
80 /** Construct object from archive_node. */
81 basic::basic(const archive_node &n, const lst &sym_lst) : flags(0), refcount(0)
83 // Reconstruct tinfo_key from class name
84 std::string class_name;
85 if (n.find_string("class", class_name))
86 tinfo_key = find_tinfo_key(class_name);
88 throw (std::runtime_error("archive node contains no class name"));
91 /** Unarchive the object. */
92 DEFAULT_UNARCHIVE(basic)
94 /** Archive the object. */
95 void basic::archive(archive_node &n) const
97 n.add_string("class", class_name());
101 // new virtual functions which can be overridden by derived classes
106 /** Output to stream.
107 * @param c print context object that describes the output formatting
108 * @param level value that is used to identify the precedence or indentation
109 * level for placing parentheses and formatting */
110 void basic::print(const print_context & c, unsigned level) const
112 if (is_of_type(c, print_tree)) {
114 c.s << std::string(level, ' ') << class_name()
115 << std::hex << ", hash=0x" << hashvalue << ", flags=0x" << flags << std::dec
116 << ", nops=" << nops()
118 for (unsigned i=0; i<nops(); ++i)
119 op(i).print(c, level + static_cast<const print_tree &>(c).delta_indent);
122 c.s << "[" << class_name() << " object]";
125 /** Little wrapper around print to be called within a debugger.
126 * This is needed because you cannot call foo.print(cout) from within the
127 * debugger because it might not know what cout is. This method can be
128 * invoked with no argument and it will simply print to stdout.
130 * @see basic::print */
131 void basic::dbgprint(void) const
133 this->print(std::cerr);
134 std::cerr << std::endl;
137 /** Little wrapper around printtree to be called within a debugger.
139 * @see basic::dbgprint
140 * @see basic::printtree */
141 void basic::dbgprinttree(void) const
143 this->print(print_tree(std::cerr));
146 /** Return relative operator precedence (for parenthizing output). */
147 unsigned basic::precedence(void) const
152 /** Create a new copy of this on the heap. One can think of this as simulating
153 * a virtual copy constructor which is needed for instance by the refcounted
154 * construction of an ex from a basic. */
155 basic * basic::duplicate() const
157 return new basic(*this);
160 /** Information about the object.
162 * @see class info_flags */
163 bool basic::info(unsigned inf) const
165 // all possible properties are false for basic objects
169 /** Number of operands/members. */
170 unsigned basic::nops() const
172 // iterating from 0 to nops() on atomic objects should be an empty loop,
173 // and accessing their elements is a range error. Container objects should
178 /** Return operand/member at position i. */
179 ex basic::op(int i) const
181 return (const_cast<basic *>(this))->let_op(i);
184 /** Return modifyable operand/member at position i. */
185 ex & basic::let_op(int i)
187 throw(std::out_of_range("op() out of range"));
190 ex basic::operator[](const ex & index) const
192 if (is_ex_exactly_of_type(index,numeric))
193 return op(ex_to<numeric>(index).to_int());
195 throw(std::invalid_argument("non-numeric indices not supported by this type"));
198 ex basic::operator[](int i) const
203 /** Test for occurrence of a pattern. An object 'has' a pattern if it matches
204 * the pattern itself or one of the children 'has' it. As a consequence
205 * (according to the definition of children) given e=x+y+z, e.has(x) is true
206 * but e.has(x+y) is false. */
207 bool basic::has(const ex & pattern) const
210 if (match(pattern, repl_lst))
212 for (unsigned i=0; i<nops(); i++)
213 if (op(i).has(pattern))
219 /** Construct new expression by applying the specified function to all
220 * sub-expressions (one level only, not recursively). */
221 ex basic::map(map_function & f) const
223 unsigned num = nops();
227 basic *copy = duplicate();
228 copy->setflag(status_flags::dynallocated);
229 copy->clearflag(status_flags::hash_calculated | status_flags::expanded);
231 for (unsigned i=0; i<num; i++)
232 e.let_op(i) = f(e.op(i));
236 /** Return degree of highest power in object s. */
237 int basic::degree(const ex & s) const
239 return is_equal(ex_to<basic>(s)) ? 1 : 0;
242 /** Return degree of lowest power in object s. */
243 int basic::ldegree(const ex & s) const
245 return is_equal(ex_to<basic>(s)) ? 1 : 0;
248 /** Return coefficient of degree n in object s. */
249 ex basic::coeff(const ex & s, int n) const
251 if (is_equal(ex_to<basic>(s)))
252 return n==1 ? _ex1 : _ex0;
254 return n==0 ? *this : _ex0;
257 /** Sort expanded expression in terms of powers of some object(s).
258 * @param s object(s) to sort in
259 * @param distributed recursive or distributed form (only used when s is a list) */
260 ex basic::collect(const ex & s, bool distributed) const
263 if (is_ex_of_type(s, lst)) {
265 // List of objects specified
269 return collect(s.op(0));
271 else if (distributed) {
273 // Get lower/upper degree of all symbols in list
278 int cnt; // current degree, 'counter'
279 ex coeff; // coefficient for degree 'cnt'
281 sym_info *si = new sym_info[num];
283 for (int i=0; i<num; i++) {
285 si[i].ldeg = si[i].cnt = this->ldegree(si[i].sym);
286 si[i].deg = this->degree(si[i].sym);
287 c = si[i].coeff = c.coeff(si[i].sym, si[i].cnt);
292 // Calculate coeff*x1^c1*...*xn^cn
294 for (int i=0; i<num; i++) {
296 y *= power(si[i].sym, cnt);
298 x += y * si[num - 1].coeff;
300 // Increment counters
304 if (si[n].cnt <= si[n].deg) {
305 // Update coefficients
311 for (int i=n; i<num; i++)
312 c = si[i].coeff = c.coeff(si[i].sym, si[i].cnt);
317 si[n].cnt = si[n].ldeg;
328 for (int n=s.nops()-1; n>=0; n--)
334 // Only one object specified
335 for (int n=this->ldegree(s); n<=this->degree(s); ++n)
336 x += this->coeff(s,n)*power(s,n);
339 // correct for lost fractional arguments and return
340 return x + (*this - x).expand();
343 /** Perform automatic non-interruptive term rewriting rules. */
344 ex basic::eval(int level) const
346 // There is nothing to do for basic objects:
350 /** Function object to be applied by basic::evalf(). */
351 struct evalf_map_function : public map_function {
353 evalf_map_function(int l) : level(l) {}
354 ex operator()(const ex & e) { return evalf(e, level); }
357 /** Evaluate object numerically. */
358 ex basic::evalf(int level) const
365 else if (level == -max_recursion_level)
366 throw(std::runtime_error("max recursion level reached"));
368 evalf_map_function map_evalf(level - 1);
369 return map(map_evalf);
374 /** Function object to be applied by basic::evalm(). */
375 struct evalm_map_function : public map_function {
376 ex operator()(const ex & e) { return evalm(e); }
379 /** Evaluate sums, products and integer powers of matrices. */
380 ex basic::evalm(void) const
385 return map(map_evalm);
388 /** Perform automatic symbolic evaluations on indexed expression that
389 * contains this object as the base expression. */
390 ex basic::eval_indexed(const basic & i) const
391 // this function can't take a "const ex & i" because that would result
392 // in an infinite eval() loop
394 // There is nothing to do for basic objects
398 /** Add two indexed expressions. They are guaranteed to be of class indexed
399 * (or a subclass) and their indices are compatible. This function is used
400 * internally by simplify_indexed().
402 * @param self First indexed expression; it's base object is *this
403 * @param other Second indexed expression
404 * @return sum of self and other
405 * @see ex::simplify_indexed() */
406 ex basic::add_indexed(const ex & self, const ex & other) const
411 /** Multiply an indexed expression with a scalar. This function is used
412 * internally by simplify_indexed().
414 * @param self Indexed expression; it's base object is *this
415 * @param other Numeric value
416 * @return product of self and other
417 * @see ex::simplify_indexed() */
418 ex basic::scalar_mul_indexed(const ex & self, const numeric & other) const
423 /** Try to contract two indexed expressions that appear in the same product.
424 * If a contraction exists, the function overwrites one or both of the
425 * expressions and returns true. Otherwise it returns false. It is
426 * guaranteed that both expressions are of class indexed (or a subclass)
427 * and that at least one dummy index has been found. This functions is
428 * used internally by simplify_indexed().
430 * @param self Pointer to first indexed expression; it's base object is *this
431 * @param other Pointer to second indexed expression
432 * @param v The complete vector of factors
433 * @return true if the contraction was successful, false otherwise
434 * @see ex::simplify_indexed() */
435 bool basic::contract_with(exvector::iterator self, exvector::iterator other, exvector & v) const
441 /** Check whether the expression matches a given pattern. For every wildcard
442 * object in the pattern, an expression of the form "wildcard == matching_expression"
443 * is added to repl_lst. */
444 bool basic::match(const ex & pattern, lst & repl_lst) const
447 Sweet sweet shapes, sweet sweet shapes,
448 That's the key thing, right right.
449 Feed feed face, feed feed shapes,
450 But who is the king tonight?
451 Who is the king tonight?
452 Pattern is the thing, the key thing-a-ling,
453 But who is the king of Pattern?
454 But who is the king, the king thing-a-ling,
455 Who is the king of Pattern?
456 Bog is the king, the king thing-a-ling,
457 Bog is the king of Pattern.
458 Ba bu-bu-bu-bu bu-bu-bu-bu-bu-bu bu-bu
459 Bog is the king of Pattern.
462 if (is_ex_exactly_of_type(pattern, wildcard)) {
464 // Wildcard matches anything, but check whether we already have found
465 // a match for that wildcard first (if so, it the earlier match must
466 // be the same expression)
467 for (unsigned i=0; i<repl_lst.nops(); i++) {
468 if (repl_lst.op(i).op(0).is_equal(pattern))
469 return is_equal(ex_to<basic>(repl_lst.op(i).op(1)));
471 repl_lst.append(pattern == *this);
476 // Expression must be of the same type as the pattern
477 if (tinfo() != ex_to<basic>(pattern).tinfo())
480 // Number of subexpressions must match
481 if (nops() != pattern.nops())
484 // No subexpressions? Then just compare the objects (there can't be
485 // wildcards in the pattern)
487 return is_equal_same_type(ex_to<basic>(pattern));
489 // Check whether attributes that are not subexpressions match
490 if (!match_same_type(ex_to<basic>(pattern)))
493 // Otherwise the subexpressions must match one-to-one
494 for (unsigned i=0; i<nops(); i++)
495 if (!op(i).match(pattern.op(i), repl_lst))
498 // Looks similar enough, match found
503 /** Substitute a set of objects by arbitrary expressions. The ex returned
504 * will already be evaluated. */
505 ex basic::subs(const lst & ls, const lst & lr, bool no_pattern) const
507 GINAC_ASSERT(ls.nops() == lr.nops());
510 for (unsigned i=0; i<ls.nops(); i++) {
511 if (is_equal(ex_to<basic>(ls.op(i))))
515 for (unsigned i=0; i<ls.nops(); i++) {
517 if (match(ex_to<basic>(ls.op(i)), repl_lst))
518 return lr.op(i).subs(repl_lst, true); // avoid infinite recursion when re-substituting the wildcards
525 /** Default interface of nth derivative ex::diff(s, n). It should be called
526 * instead of ::derivative(s) for first derivatives and for nth derivatives it
527 * just recurses down.
529 * @param s symbol to differentiate in
530 * @param nth order of differentiation
532 ex basic::diff(const symbol & s, unsigned nth) const
534 // trivial: zeroth derivative
538 // evaluate unevaluated *this before differentiating
539 if (!(flags & status_flags::evaluated))
540 return ex(*this).diff(s, nth);
542 ex ndiff = this->derivative(s);
543 while (!ndiff.is_zero() && // stop differentiating zeros
545 ndiff = ndiff.diff(s);
551 /** Return a vector containing the free indices of an expression. */
552 exvector basic::get_free_indices(void) const
554 return exvector(); // return an empty exvector
557 ex basic::simplify_ncmul(const exvector & v) const
559 return simplified_ncmul(v);
564 /** Function object to be applied by basic::derivative(). */
565 struct derivative_map_function : public map_function {
567 derivative_map_function(const symbol &sym) : s(sym) {}
568 ex operator()(const ex & e) { return diff(e, s); }
571 /** Default implementation of ex::diff(). It maps the operation on the
572 * operands (or returns 0 when the object has no operands).
575 ex basic::derivative(const symbol & s) const
580 derivative_map_function map_derivative(s);
581 return map(map_derivative);
585 /** Returns order relation between two objects of same type. This needs to be
586 * implemented by each class. It may never return anything else than 0,
587 * signalling equality, or +1 and -1 signalling inequality and determining
588 * the canonical ordering. (Perl hackers will wonder why C++ doesn't feature
589 * the spaceship operator <=> for denoting just this.) */
590 int basic::compare_same_type(const basic & other) const
592 return compare_pointers(this, &other);
595 /** Returns true if two objects of same type are equal. Normally needs
596 * not be reimplemented as long as it wasn't overwritten by some parent
597 * class, since it just calls compare_same_type(). The reason why this
598 * function exists is that sometimes it is easier to determine equality
599 * than an order relation and then it can be overridden. */
600 bool basic::is_equal_same_type(const basic & other) const
602 return compare_same_type(other)==0;
605 /** Returns true if the attributes of two objects are similar enough for
606 * a match. This function must not match subexpressions (this is already
607 * done by basic::match()). Only attributes not accessible by op() should
608 * be compared. This is also the reason why this function doesn't take the
609 * wildcard replacement list from match() as an argument: only subexpressions
610 * are subject to wildcard matches. Also, this function only needs to be
611 * implemented for container classes because is_equal_same_type() is
612 * automatically used instead of match_same_type() if nops() == 0.
614 * @see basic::match */
615 bool basic::match_same_type(const basic & other) const
617 // The default is to only consider subexpressions, but not any other
622 unsigned basic::return_type(void) const
624 return return_types::commutative;
627 unsigned basic::return_type_tinfo(void) const
632 /** Compute the hash value of an object and if it makes sense to store it in
633 * the objects status_flags, do so. The method inherited from class basic
634 * computes a hash value based on the type and hash values of possible
635 * members. For this reason it is well suited for container classes but
636 * atomic classes should override this implementation because otherwise they
637 * would all end up with the same hashvalue. */
638 unsigned basic::calchash(void) const
640 unsigned v = golden_ratio_hash(tinfo());
641 for (unsigned i=0; i<nops(); i++) {
642 v = rotate_left_31(v);
643 v ^= (const_cast<basic *>(this))->op(i).gethash();
646 // mask out numeric hashes:
649 // store calculated hash value only if object is already evaluated
650 if (flags & status_flags::evaluated) {
651 setflag(status_flags::hash_calculated);
658 /** Function object to be applied by basic::expand(). */
659 struct expand_map_function : public map_function {
661 expand_map_function(unsigned o) : options(o) {}
662 ex operator()(const ex & e) { return expand(e, options); }
665 /** Expand expression, i.e. multiply it out and return the result as a new
667 ex basic::expand(unsigned options) const
670 return (options == 0) ? setflag(status_flags::expanded) : *this;
672 expand_map_function map_expand(options);
673 return ex_to<basic>(map(map_expand)).setflag(options == 0 ? status_flags::expanded : 0);
679 // non-virtual functions in this class
684 /** Substitute objects in an expression (syntactic substitution) and return
685 * the result as a new expression. There are two valid types of
686 * replacement arguments: 1) a relational like object==ex and 2) a list of
687 * relationals lst(object1==ex1,object2==ex2,...), which is converted to
688 * subs(lst(object1,object2,...),lst(ex1,ex2,...)). */
689 ex basic::subs(const ex & e, bool no_pattern) const
691 if (e.info(info_flags::relation_equal)) {
692 return subs(lst(e), no_pattern);
694 if (!e.info(info_flags::list)) {
695 throw(std::invalid_argument("basic::subs(ex): argument must be a list"));
699 for (unsigned i=0; i<e.nops(); i++) {
701 if (!r.info(info_flags::relation_equal)) {
702 throw(std::invalid_argument("basic::subs(ex): argument must be a list of equations"));
707 return subs(ls, lr, no_pattern);
710 /** Compare objects to establish canonical ordering.
711 * All compare functions return: -1 for *this less than other, 0 equal,
713 int basic::compare(const basic & other) const
715 unsigned hash_this = gethash();
716 unsigned hash_other = other.gethash();
718 if (hash_this<hash_other) return -1;
719 if (hash_this>hash_other) return 1;
721 unsigned typeid_this = tinfo();
722 unsigned typeid_other = other.tinfo();
724 if (typeid_this<typeid_other) {
725 // std::cout << "hash collision, different types: "
726 // << *this << " and " << other << std::endl;
727 // this->print(print_tree(std::cout));
728 // std::cout << " and ";
729 // other.print(print_tree(std::cout));
730 // std::cout << std::endl;
733 if (typeid_this>typeid_other) {
734 // std::cout << "hash collision, different types: "
735 // << *this << " and " << other << std::endl;
736 // this->print(print_tree(std::cout));
737 // std::cout << " and ";
738 // other.print(print_tree(std::cout));
739 // std::cout << std::endl;
743 GINAC_ASSERT(typeid(*this)==typeid(other));
745 // int cmpval = compare_same_type(other);
746 // if ((cmpval!=0) && (hash_this<0x80000000U)) {
747 // std::cout << "hash collision, same type: "
748 // << *this << " and " << other << std::endl;
749 // this->print(print_tree(std::cout));
750 // std::cout << " and ";
751 // other.print(print_tree(std::cout));
752 // std::cout << std::endl;
756 return compare_same_type(other);
759 /** Test for equality.
760 * This is only a quick test, meaning objects should be in the same domain.
761 * You might have to .expand(), .normal() objects first, depending on the
762 * domain of your computation, to get a more reliable answer.
764 * @see is_equal_same_type */
765 bool basic::is_equal(const basic & other) const
767 if (this->gethash()!=other.gethash())
769 if (this->tinfo()!=other.tinfo())
772 GINAC_ASSERT(typeid(*this)==typeid(other));
774 return is_equal_same_type(other);
779 /** Stop further evaluation.
781 * @see basic::eval */
782 const basic & basic::hold(void) const
784 return setflag(status_flags::evaluated);
787 /** Ensure the object may be modified without hurting others, throws if this
788 * is not the case. */
789 void basic::ensure_if_modifiable(void) const
791 if (this->refcount>1)
792 throw(std::runtime_error("cannot modify multiply referenced object"));
793 clearflag(status_flags::hash_calculated);
800 int max_recursion_level = 1024;