3 * Implementation of GiNaC's ABC. */
6 * GiNaC Copyright (C) 1999-2001 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
242 /** Return degree of lowest power in object s. */
243 int basic::ldegree(const ex & s) const
248 /** Return coefficient of degree n in object s. */
249 ex basic::coeff(const ex & s, int n) const
251 return n==0 ? *this : _ex0;
254 /** Sort expanded expression in terms of powers of some object(s).
255 * @param s object(s) to sort in
256 * @param distributed recursive or distributed form (only used when s is a list) */
257 ex basic::collect(const ex & s, bool distributed) const
260 if (is_ex_of_type(s, lst)) {
262 // List of objects specified
266 return collect(s.op(0));
268 else if (distributed) {
270 // Get lower/upper degree of all symbols in list
275 int cnt; // current degree, 'counter'
276 ex coeff; // coefficient for degree 'cnt'
278 sym_info *si = new sym_info[num];
280 for (int i=0; i<num; i++) {
282 si[i].ldeg = si[i].cnt = this->ldegree(si[i].sym);
283 si[i].deg = this->degree(si[i].sym);
284 c = si[i].coeff = c.coeff(si[i].sym, si[i].cnt);
289 // Calculate coeff*x1^c1*...*xn^cn
291 for (int i=0; i<num; i++) {
293 y *= power(si[i].sym, cnt);
295 x += y * si[num - 1].coeff;
297 // Increment counters
301 if (si[n].cnt <= si[n].deg) {
302 // Update coefficients
308 for (int i=n; i<num; i++)
309 c = si[i].coeff = c.coeff(si[i].sym, si[i].cnt);
314 si[n].cnt = si[n].ldeg;
325 for (int n=s.nops()-1; n>=0; n--)
331 // Only one object specified
332 for (int n=this->ldegree(s); n<=this->degree(s); ++n)
333 x += this->coeff(s,n)*power(s,n);
336 // correct for lost fractional arguments and return
337 return x + (*this - x).expand();
340 /** Perform automatic non-interruptive term rewriting rules. */
341 ex basic::eval(int level) const
343 // There is nothing to do for basic objects:
347 /** Function object to be applied by basic::evalf(). */
348 struct evalf_map_function : public map_function {
350 evalf_map_function(int l) : level(l) {}
351 ex operator()(const ex & e) { return evalf(e, level); }
354 /** Evaluate object numerically. */
355 ex basic::evalf(int level) const
362 else if (level == -max_recursion_level)
363 throw(std::runtime_error("max recursion level reached"));
365 evalf_map_function map_evalf(level - 1);
366 return map(map_evalf);
371 /** Function object to be applied by basic::evalm(). */
372 struct evalm_map_function : public map_function {
373 ex operator()(const ex & e) { return evalm(e); }
376 /** Evaluate sums, products and integer powers of matrices. */
377 ex basic::evalm(void) const
382 return map(map_evalm);
385 /** Perform automatic symbolic evaluations on indexed expression that
386 * contains this object as the base expression. */
387 ex basic::eval_indexed(const basic & i) const
388 // this function can't take a "const ex & i" because that would result
389 // in an infinite eval() loop
391 // There is nothing to do for basic objects
395 /** Add two indexed expressions. They are guaranteed to be of class indexed
396 * (or a subclass) and their indices are compatible. This function is used
397 * internally by simplify_indexed().
399 * @param self First indexed expression; it's base object is *this
400 * @param other Second indexed expression
401 * @return sum of self and other
402 * @see ex::simplify_indexed() */
403 ex basic::add_indexed(const ex & self, const ex & other) const
408 /** Multiply an indexed expression with a scalar. This function is used
409 * internally by simplify_indexed().
411 * @param self Indexed expression; it's base object is *this
412 * @param other Numeric value
413 * @return product of self and other
414 * @see ex::simplify_indexed() */
415 ex basic::scalar_mul_indexed(const ex & self, const numeric & other) const
420 /** Try to contract two indexed expressions that appear in the same product.
421 * If a contraction exists, the function overwrites one or both of the
422 * expressions and returns true. Otherwise it returns false. It is
423 * guaranteed that both expressions are of class indexed (or a subclass)
424 * and that at least one dummy index has been found. This functions is
425 * used internally by simplify_indexed().
427 * @param self Pointer to first indexed expression; it's base object is *this
428 * @param other Pointer to second indexed expression
429 * @param v The complete vector of factors
430 * @return true if the contraction was successful, false otherwise
431 * @see ex::simplify_indexed() */
432 bool basic::contract_with(exvector::iterator self, exvector::iterator other, exvector & v) const
438 /** Check whether the expression matches a given pattern. For every wildcard
439 * object in the pattern, an expression of the form "wildcard == matching_expression"
440 * is added to repl_lst. */
441 bool basic::match(const ex & pattern, lst & repl_lst) const
444 Sweet sweet shapes, sweet sweet shapes,
445 That's the key thing, right right.
446 Feed feed face, feed feed shapes,
447 But who is the king tonight?
448 Who is the king tonight?
449 Pattern is the thing, the key thing-a-ling,
450 But who is the king of Pattern?
451 But who is the king, the king thing-a-ling,
452 Who is the king of Pattern?
453 Bog is the king, the king thing-a-ling,
454 Bog is the king of Pattern.
455 Ba bu-bu-bu-bu bu-bu-bu-bu-bu-bu bu-bu
456 Bog is the king of Pattern.
459 if (is_ex_exactly_of_type(pattern, wildcard)) {
461 // Wildcard matches anything, but check whether we already have found
462 // a match for that wildcard first (if so, it the earlier match must
463 // be the same expression)
464 for (unsigned i=0; i<repl_lst.nops(); i++) {
465 if (repl_lst.op(i).op(0).is_equal(pattern))
466 return is_equal(ex_to<basic>(repl_lst.op(i).op(1)));
468 repl_lst.append(pattern == *this);
473 // Expression must be of the same type as the pattern
474 if (tinfo() != ex_to<basic>(pattern).tinfo())
477 // Number of subexpressions must match
478 if (nops() != pattern.nops())
481 // No subexpressions? Then just compare the objects (there can't be
482 // wildcards in the pattern)
484 return is_equal_same_type(ex_to<basic>(pattern));
486 // Check whether attributes that are not subexpressions match
487 if (!match_same_type(ex_to<basic>(pattern)))
490 // Otherwise the subexpressions must match one-to-one
491 for (unsigned i=0; i<nops(); i++)
492 if (!op(i).match(pattern.op(i), repl_lst))
495 // Looks similar enough, match found
500 /** Substitute a set of objects by arbitrary expressions. The ex returned
501 * will already be evaluated. */
502 ex basic::subs(const lst & ls, const lst & lr, bool no_pattern) const
504 GINAC_ASSERT(ls.nops() == lr.nops());
507 for (unsigned i=0; i<ls.nops(); i++) {
508 if (is_equal(ex_to<basic>(ls.op(i))))
512 for (unsigned i=0; i<ls.nops(); i++) {
514 if (match(ex_to<basic>(ls.op(i)), repl_lst))
515 return lr.op(i).subs(repl_lst, true); // avoid infinite recursion when re-substituting the wildcards
522 /** Default interface of nth derivative ex::diff(s, n). It should be called
523 * instead of ::derivative(s) for first derivatives and for nth derivatives it
524 * just recurses down.
526 * @param s symbol to differentiate in
527 * @param nth order of differentiation
529 ex basic::diff(const symbol & s, unsigned nth) const
531 // trivial: zeroth derivative
535 // evaluate unevaluated *this before differentiating
536 if (!(flags & status_flags::evaluated))
537 return ex(*this).diff(s, nth);
539 ex ndiff = this->derivative(s);
540 while (!ndiff.is_zero() && // stop differentiating zeros
542 ndiff = ndiff.diff(s);
548 /** Return a vector containing the free indices of an expression. */
549 exvector basic::get_free_indices(void) const
551 return exvector(); // return an empty exvector
554 ex basic::simplify_ncmul(const exvector & v) const
556 return simplified_ncmul(v);
561 /** Function object to be applied by basic::derivative(). */
562 struct derivative_map_function : public map_function {
564 derivative_map_function(const symbol &sym) : s(sym) {}
565 ex operator()(const ex & e) { return diff(e, s); }
568 /** Default implementation of ex::diff(). It maps the operation on the
569 * operands (or returns 0 when the object has no operands).
572 ex basic::derivative(const symbol & s) const
577 derivative_map_function map_derivative(s);
578 return map(map_derivative);
582 /** Returns order relation between two objects of same type. This needs to be
583 * implemented by each class. It may never return anything else than 0,
584 * signalling equality, or +1 and -1 signalling inequality and determining
585 * the canonical ordering. (Perl hackers will wonder why C++ doesn't feature
586 * the spaceship operator <=> for denoting just this.) */
587 int basic::compare_same_type(const basic & other) const
589 return compare_pointers(this, &other);
592 /** Returns true if two objects of same type are equal. Normally needs
593 * not be reimplemented as long as it wasn't overwritten by some parent
594 * class, since it just calls compare_same_type(). The reason why this
595 * function exists is that sometimes it is easier to determine equality
596 * than an order relation and then it can be overridden. */
597 bool basic::is_equal_same_type(const basic & other) const
599 return compare_same_type(other)==0;
602 /** Returns true if the attributes of two objects are similar enough for
603 * a match. This function must not match subexpressions (this is already
604 * done by basic::match()). Only attributes not accessible by op() should
605 * be compared. This is also the reason why this function doesn't take the
606 * wildcard replacement list from match() as an argument: only subexpressions
607 * are subject to wildcard matches. Also, this function only needs to be
608 * implemented for container classes because is_equal_same_type() is
609 * automatically used instead of match_same_type() if nops() == 0.
611 * @see basic::match */
612 bool basic::match_same_type(const basic & other) const
614 // The default is to only consider subexpressions, but not any other
619 unsigned basic::return_type(void) const
621 return return_types::commutative;
624 unsigned basic::return_type_tinfo(void) const
629 /** Compute the hash value of an object and if it makes sense to store it in
630 * the objects status_flags, do so. The method inherited from class basic
631 * computes a hash value based on the type and hash values of possible
632 * members. For this reason it is well suited for container classes but
633 * atomic classes should override this implementation because otherwise they
634 * would all end up with the same hashvalue. */
635 unsigned basic::calchash(void) const
637 unsigned v = golden_ratio_hash(tinfo());
638 for (unsigned i=0; i<nops(); i++) {
639 v = rotate_left_31(v);
640 v ^= (const_cast<basic *>(this))->op(i).gethash();
643 // mask out numeric hashes:
646 // store calculated hash value only if object is already evaluated
647 if (flags & status_flags::evaluated) {
648 setflag(status_flags::hash_calculated);
655 /** Function object to be applied by basic::expand(). */
656 struct expand_map_function : public map_function {
658 expand_map_function(unsigned o) : options(o) {}
659 ex operator()(const ex & e) { return expand(e, options); }
662 /** Expand expression, i.e. multiply it out and return the result as a new
664 ex basic::expand(unsigned options) const
667 return (options == 0) ? setflag(status_flags::expanded) : *this;
669 expand_map_function map_expand(options);
670 return ex_to<basic>(map(map_expand)).setflag(options == 0 ? status_flags::expanded : 0);
676 // non-virtual functions in this class
681 /** Substitute objects in an expression (syntactic substitution) and return
682 * the result as a new expression. There are two valid types of
683 * replacement arguments: 1) a relational like object==ex and 2) a list of
684 * relationals lst(object1==ex1,object2==ex2,...), which is converted to
685 * subs(lst(object1,object2,...),lst(ex1,ex2,...)). */
686 ex basic::subs(const ex & e, bool no_pattern) const
688 if (e.info(info_flags::relation_equal)) {
689 return subs(lst(e), no_pattern);
691 if (!e.info(info_flags::list)) {
692 throw(std::invalid_argument("basic::subs(ex): argument must be a list"));
696 for (unsigned i=0; i<e.nops(); i++) {
698 if (!r.info(info_flags::relation_equal)) {
699 throw(std::invalid_argument("basic::subs(ex): argument must be a list of equations"));
704 return subs(ls, lr, no_pattern);
707 /** Compare objects to establish canonical ordering.
708 * All compare functions return: -1 for *this less than other, 0 equal,
710 int basic::compare(const basic & other) const
712 unsigned hash_this = gethash();
713 unsigned hash_other = other.gethash();
715 if (hash_this<hash_other) return -1;
716 if (hash_this>hash_other) return 1;
718 unsigned typeid_this = tinfo();
719 unsigned typeid_other = other.tinfo();
721 if (typeid_this<typeid_other) {
722 // std::cout << "hash collision, different types: "
723 // << *this << " and " << other << std::endl;
724 // this->print(print_tree(std::cout));
725 // std::cout << " and ";
726 // other.print(print_tree(std::cout));
727 // std::cout << std::endl;
730 if (typeid_this>typeid_other) {
731 // std::cout << "hash collision, different types: "
732 // << *this << " and " << other << std::endl;
733 // this->print(print_tree(std::cout));
734 // std::cout << " and ";
735 // other.print(print_tree(std::cout));
736 // std::cout << std::endl;
740 GINAC_ASSERT(typeid(*this)==typeid(other));
742 // int cmpval = compare_same_type(other);
743 // if ((cmpval!=0) && (hash_this<0x80000000U)) {
744 // std::cout << "hash collision, same type: "
745 // << *this << " and " << other << std::endl;
746 // this->print(print_tree(std::cout));
747 // std::cout << " and ";
748 // other.print(print_tree(std::cout));
749 // std::cout << std::endl;
753 return compare_same_type(other);
756 /** Test for equality.
757 * This is only a quick test, meaning objects should be in the same domain.
758 * You might have to .expand(), .normal() objects first, depending on the
759 * domain of your computation, to get a more reliable answer.
761 * @see is_equal_same_type */
762 bool basic::is_equal(const basic & other) const
764 if (this->gethash()!=other.gethash())
766 if (this->tinfo()!=other.tinfo())
769 GINAC_ASSERT(typeid(*this)==typeid(other));
771 return is_equal_same_type(other);
776 /** Stop further evaluation.
778 * @see basic::eval */
779 const basic & basic::hold(void) const
781 return setflag(status_flags::evaluated);
784 /** Ensure the object may be modified without hurting others, throws if this
785 * is not the case. */
786 void basic::ensure_if_modifiable(void) const
788 if (this->refcount>1)
789 throw(std::runtime_error("cannot modify multiply referenced object"));
790 clearflag(status_flags::hash_calculated);
797 int max_recursion_level = 1024;