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"
45 GINAC_IMPLEMENT_REGISTERED_CLASS_NO_CTORS(basic, void)
48 // default ctor, dtor, copy ctor assignment operator and helpers
53 basic::basic(const basic & other) : tinfo_key(TINFO_basic), flags(0), refcount(0)
55 debugmsg("basic copy ctor", LOGLEVEL_CONSTRUCT);
59 const basic & basic::operator=(const basic & other)
61 debugmsg("basic operator=", LOGLEVEL_ASSIGNMENT);
71 // none (all conditionally inlined)
77 // none (all conditionally inlined)
83 /** Construct object from archive_node. */
84 basic::basic(const archive_node &n, const lst &sym_lst) : flags(0), refcount(0)
86 debugmsg("basic ctor from archive_node", LOGLEVEL_CONSTRUCT);
88 // Reconstruct tinfo_key from class name
89 std::string class_name;
90 if (n.find_string("class", class_name))
91 tinfo_key = find_tinfo_key(class_name);
93 throw (std::runtime_error("archive node contains no class name"));
96 /** Unarchive the object. */
97 DEFAULT_UNARCHIVE(basic)
99 /** Archive the object. */
100 void basic::archive(archive_node &n) const
102 n.add_string("class", class_name());
106 // new virtual functions which can be overridden by derived classes
111 /** Output to stream.
112 * @param c print context object that describes the output formatting
113 * @param level value that is used to identify the precedence or indentation
114 * level for placing parentheses and formatting */
115 void basic::print(const print_context & c, unsigned level) const
117 debugmsg("basic print", LOGLEVEL_PRINT);
119 if (is_of_type(c, print_tree)) {
121 c.s << std::string(level, ' ') << class_name()
122 << std::hex << ", hash=0x" << hashvalue << ", flags=0x" << flags << std::dec
123 << ", nops=" << nops()
125 for (unsigned i=0; i<nops(); ++i)
126 op(i).print(c, level + static_cast<const print_tree &>(c).delta_indent);
129 c.s << "[" << class_name() << " object]";
132 /** Little wrapper around print to be called within a debugger.
133 * This is needed because you cannot call foo.print(cout) from within the
134 * debugger because it might not know what cout is. This method can be
135 * invoked with no argument and it will simply print to stdout.
137 * @see basic::print */
138 void basic::dbgprint(void) const
140 this->print(std::cerr);
141 std::cerr << std::endl;
144 /** Little wrapper around printtree to be called within a debugger.
146 * @see basic::dbgprint
147 * @see basic::printtree */
148 void basic::dbgprinttree(void) const
150 this->print(print_tree(std::cerr));
153 /** Return relative operator precedence (for parenthizing output). */
154 unsigned basic::precedence(void) const
159 /** Create a new copy of this on the heap. One can think of this as simulating
160 * a virtual copy constructor which is needed for instance by the refcounted
161 * construction of an ex from a basic. */
162 basic * basic::duplicate() const
164 debugmsg("basic duplicate",LOGLEVEL_DUPLICATE);
165 return new basic(*this);
168 /** Information about the object.
170 * @see class info_flags */
171 bool basic::info(unsigned inf) const
173 // all possible properties are false for basic objects
177 /** Number of operands/members. */
178 unsigned basic::nops() const
180 // iterating from 0 to nops() on atomic objects should be an empty loop,
181 // and accessing their elements is a range error. Container objects should
186 /** Return operand/member at position i. */
187 ex basic::op(int i) const
189 return (const_cast<basic *>(this))->let_op(i);
192 /** Return modifyable operand/member at position i. */
193 ex & basic::let_op(int i)
195 throw(std::out_of_range("op() out of range"));
198 ex basic::operator[](const ex & index) const
200 if (is_ex_exactly_of_type(index,numeric))
201 return op(ex_to<numeric>(index).to_int());
203 throw(std::invalid_argument("non-numeric indices not supported by this type"));
206 ex basic::operator[](int i) const
211 /** Test for occurrence of a pattern. An object 'has' a pattern if it matches
212 * the pattern itself or one of the children 'has' it. As a consequence
213 * (according to the definition of children) given e=x+y+z, e.has(x) is true
214 * but e.has(x+y) is false. */
215 bool basic::has(const ex & pattern) const
218 if (match(pattern, repl_lst))
220 for (unsigned i=0; i<nops(); i++)
221 if (op(i).has(pattern))
227 /** Construct new expression by applying the specified function to all
228 * sub-expressions (one level only, not recursively). */
229 ex basic::map(map_function & f) const
231 unsigned num = nops();
235 basic *copy = duplicate();
236 copy->setflag(status_flags::dynallocated);
237 copy->clearflag(status_flags::hash_calculated | status_flags::expanded);
239 for (unsigned i=0; i<num; i++)
240 e.let_op(i) = f(e.op(i));
244 /** Return degree of highest power in object s. */
245 int basic::degree(const ex & s) const
250 /** Return degree of lowest power in object s. */
251 int basic::ldegree(const ex & s) const
256 /** Return coefficient of degree n in object s. */
257 ex basic::coeff(const ex & s, int n) const
259 return n==0 ? *this : _ex0();
262 /** Sort expanded expression in terms of powers of some object(s).
263 * @param s object(s) to sort in
264 * @param distributed recursive or distributed form (only used when s is a list) */
265 ex basic::collect(const ex & s, bool distributed) const
268 if (is_ex_of_type(s, lst)) {
270 // List of objects specified
274 return collect(s.op(0));
276 else if (distributed) {
278 // Get lower/upper degree of all symbols in list
283 int cnt; // current degree, 'counter'
284 ex coeff; // coefficient for degree 'cnt'
286 sym_info *si = new sym_info[num];
288 for (int i=0; i<num; i++) {
290 si[i].ldeg = si[i].cnt = this->ldegree(si[i].sym);
291 si[i].deg = this->degree(si[i].sym);
292 c = si[i].coeff = c.coeff(si[i].sym, si[i].cnt);
297 // Calculate coeff*x1^c1*...*xn^cn
299 for (int i=0; i<num; i++) {
301 y *= power(si[i].sym, cnt);
303 x += y * si[num - 1].coeff;
305 // Increment counters
309 if (si[n].cnt <= si[n].deg) {
310 // Update coefficients
316 for (int i=n; i<num; i++)
317 c = si[i].coeff = c.coeff(si[i].sym, si[i].cnt);
322 si[n].cnt = si[n].ldeg;
333 for (int n=s.nops()-1; n>=0; n--)
339 // Only one object specified
340 for (int n=this->ldegree(s); n<=this->degree(s); ++n)
341 x += this->coeff(s,n)*power(s,n);
344 // correct for lost fractional arguments and return
345 return x + (*this - x).expand();
348 /** Perform automatic non-interruptive term rewriting rules. */
349 ex basic::eval(int level) const
351 // There is nothing to do for basic objects:
355 /** Function object to be applied by basic::evalf(). */
356 struct evalf_map_function : public map_function {
358 evalf_map_function(int l) : level(l) {}
359 ex operator()(const ex & e) { return evalf(e, level); }
362 /** Evaluate object numerically. */
363 ex basic::evalf(int level) const
370 else if (level == -max_recursion_level)
371 throw(std::runtime_error("max recursion level reached"));
373 evalf_map_function map_evalf(level - 1);
374 return map(map_evalf);
379 /** Function object to be applied by basic::evalm(). */
380 struct evalm_map_function : public map_function {
381 ex operator()(const ex & e) { return evalm(e); }
384 /** Evaluate sums, products and integer powers of matrices. */
385 ex basic::evalm(void) const
390 return map(map_evalm);
393 /** Perform automatic symbolic evaluations on indexed expression that
394 * contains this object as the base expression. */
395 ex basic::eval_indexed(const basic & i) const
396 // this function can't take a "const ex & i" because that would result
397 // in an infinite eval() loop
399 // There is nothing to do for basic objects
403 /** Add two indexed expressions. They are guaranteed to be of class indexed
404 * (or a subclass) and their indices are compatible. This function is used
405 * internally by simplify_indexed().
407 * @param self First indexed expression; it's base object is *this
408 * @param other Second indexed expression
409 * @return sum of self and other
410 * @see ex::simplify_indexed() */
411 ex basic::add_indexed(const ex & self, const ex & other) const
416 /** Multiply an indexed expression with a scalar. This function is used
417 * internally by simplify_indexed().
419 * @param self Indexed expression; it's base object is *this
420 * @param other Numeric value
421 * @return product of self and other
422 * @see ex::simplify_indexed() */
423 ex basic::scalar_mul_indexed(const ex & self, const numeric & other) const
428 /** Try to contract two indexed expressions that appear in the same product.
429 * If a contraction exists, the function overwrites one or both of the
430 * expressions and returns true. Otherwise it returns false. It is
431 * guaranteed that both expressions are of class indexed (or a subclass)
432 * and that at least one dummy index has been found. This functions is
433 * used internally by simplify_indexed().
435 * @param self Pointer to first indexed expression; it's base object is *this
436 * @param other Pointer to second indexed expression
437 * @param v The complete vector of factors
438 * @return true if the contraction was successful, false otherwise
439 * @see ex::simplify_indexed() */
440 bool basic::contract_with(exvector::iterator self, exvector::iterator other, exvector & v) const
446 /** Check whether the expression matches a given pattern. For every wildcard
447 * object in the pattern, an expression of the form "wildcard == matching_expression"
448 * is added to repl_lst. */
449 bool basic::match(const ex & pattern, lst & repl_lst) const
452 Sweet sweet shapes, sweet sweet shapes,
453 That's the key thing, right right.
454 Feed feed face, feed feed shapes,
455 But who is the king tonight?
456 Who is the king tonight?
457 Pattern is the thing, the key thing-a-ling,
458 But who is the king of Pattern?
459 But who is the king, the king thing-a-ling,
460 Who is the king of Pattern?
461 Bog is the king, the king thing-a-ling,
462 Bog is the king of Pattern.
463 Ba bu-bu-bu-bu bu-bu-bu-bu-bu-bu bu-bu
464 Bog is the king of Pattern.
467 if (is_ex_exactly_of_type(pattern, wildcard)) {
469 // Wildcard matches anything, but check whether we already have found
470 // a match for that wildcard first (if so, it the earlier match must
471 // be the same expression)
472 for (unsigned i=0; i<repl_lst.nops(); i++) {
473 if (repl_lst.op(i).op(0).is_equal(pattern))
474 return is_equal(ex_to<basic>(repl_lst.op(i).op(1)));
476 repl_lst.append(pattern == *this);
481 // Expression must be of the same type as the pattern
482 if (tinfo() != ex_to<basic>(pattern).tinfo())
485 // Number of subexpressions must match
486 if (nops() != pattern.nops())
489 // No subexpressions? Then just compare the objects (there can't be
490 // wildcards in the pattern)
492 return is_equal_same_type(ex_to<basic>(pattern));
494 // Check whether attributes that are not subexpressions match
495 if (!match_same_type(ex_to<basic>(pattern)))
498 // Otherwise the subexpressions must match one-to-one
499 for (unsigned i=0; i<nops(); i++)
500 if (!op(i).match(pattern.op(i), repl_lst))
503 // Looks similar enough, match found
508 /** Substitute a set of objects by arbitrary expressions. The ex returned
509 * will already be evaluated. */
510 ex basic::subs(const lst & ls, const lst & lr, bool no_pattern) const
512 GINAC_ASSERT(ls.nops() == lr.nops());
515 for (unsigned i=0; i<ls.nops(); i++) {
516 if (is_equal(ex_to<basic>(ls.op(i))))
520 for (unsigned i=0; i<ls.nops(); i++) {
522 if (match(ex_to<basic>(ls.op(i)), repl_lst))
523 return lr.op(i).subs(repl_lst, true); // avoid infinite recursion when re-substituting the wildcards
530 /** Default interface of nth derivative ex::diff(s, n). It should be called
531 * instead of ::derivative(s) for first derivatives and for nth derivatives it
532 * just recurses down.
534 * @param s symbol to differentiate in
535 * @param nth order of differentiation
537 ex basic::diff(const symbol & s, unsigned nth) const
539 // trivial: zeroth derivative
543 // evaluate unevaluated *this before differentiating
544 if (!(flags & status_flags::evaluated))
545 return ex(*this).diff(s, nth);
547 ex ndiff = this->derivative(s);
548 while (!ndiff.is_zero() && // stop differentiating zeros
550 ndiff = ndiff.diff(s);
556 /** Return a vector containing the free indices of an expression. */
557 exvector basic::get_free_indices(void) const
559 return exvector(); // return an empty exvector
562 ex basic::simplify_ncmul(const exvector & v) const
564 return simplified_ncmul(v);
569 /** Function object to be applied by basic::derivative(). */
570 struct derivative_map_function : public map_function {
572 derivative_map_function(const symbol &sym) : s(sym) {}
573 ex operator()(const ex & e) { return diff(e, s); }
576 /** Default implementation of ex::diff(). It maps the operation on the
577 * operands (or returns 0 when the object has no operands).
580 ex basic::derivative(const symbol & s) const
585 derivative_map_function map_derivative(s);
586 return map(map_derivative);
590 /** Returns order relation between two objects of same type. This needs to be
591 * implemented by each class. It may never return anything else than 0,
592 * signalling equality, or +1 and -1 signalling inequality and determining
593 * the canonical ordering. (Perl hackers will wonder why C++ doesn't feature
594 * the spaceship operator <=> for denoting just this.) */
595 int basic::compare_same_type(const basic & other) const
597 return compare_pointers(this, &other);
600 /** Returns true if two objects of same type are equal. Normally needs
601 * not be reimplemented as long as it wasn't overwritten by some parent
602 * class, since it just calls compare_same_type(). The reason why this
603 * function exists is that sometimes it is easier to determine equality
604 * than an order relation and then it can be overridden. */
605 bool basic::is_equal_same_type(const basic & other) const
607 return compare_same_type(other)==0;
610 /** Returns true if the attributes of two objects are similar enough for
611 * a match. This function must not match subexpressions (this is already
612 * done by basic::match()). Only attributes not accessible by op() should
613 * be compared. This is also the reason why this function doesn't take the
614 * wildcard replacement list from match() as an argument: only subexpressions
615 * are subject to wildcard matches. Also, this function only needs to be
616 * implemented for container classes because is_equal_same_type() is
617 * automatically used instead of match_same_type() if nops() == 0.
619 * @see basic::match */
620 bool basic::match_same_type(const basic & other) const
622 // The default is to only consider subexpressions, but not any other
627 unsigned basic::return_type(void) const
629 return return_types::commutative;
632 unsigned basic::return_type_tinfo(void) const
637 /** Compute the hash value of an object and if it makes sense to store it in
638 * the objects status_flags, do so. The method inherited from class basic
639 * computes a hash value based on the type and hash values of possible
640 * members. For this reason it is well suited for container classes but
641 * atomic classes should override this implementation because otherwise they
642 * would all end up with the same hashvalue. */
643 unsigned basic::calchash(void) const
645 unsigned v = golden_ratio_hash(tinfo());
646 for (unsigned i=0; i<nops(); i++) {
647 v = rotate_left_31(v);
648 v ^= (const_cast<basic *>(this))->op(i).gethash();
651 // mask out numeric hashes:
654 // store calculated hash value only if object is already evaluated
655 if (flags & status_flags::evaluated) {
656 setflag(status_flags::hash_calculated);
663 /** Function object to be applied by basic::expand(). */
664 struct expand_map_function : public map_function {
666 expand_map_function(unsigned o) : options(o) {}
667 ex operator()(const ex & e) { return expand(e, options); }
670 /** Expand expression, i.e. multiply it out and return the result as a new
672 ex basic::expand(unsigned options) const
675 return (options == 0) ? setflag(status_flags::expanded) : *this;
677 expand_map_function map_expand(options);
678 return ex_to<basic>(map(map_expand)).setflag(options == 0 ? status_flags::expanded : 0);
684 // non-virtual functions in this class
689 /** Substitute objects in an expression (syntactic substitution) and return
690 * the result as a new expression. There are two valid types of
691 * replacement arguments: 1) a relational like object==ex and 2) a list of
692 * relationals lst(object1==ex1,object2==ex2,...), which is converted to
693 * subs(lst(object1,object2,...),lst(ex1,ex2,...)). */
694 ex basic::subs(const ex & e, bool no_pattern) const
696 if (e.info(info_flags::relation_equal)) {
697 return subs(lst(e), no_pattern);
699 if (!e.info(info_flags::list)) {
700 throw(std::invalid_argument("basic::subs(ex): argument must be a list"));
704 for (unsigned i=0; i<e.nops(); i++) {
706 if (!r.info(info_flags::relation_equal)) {
707 throw(std::invalid_argument("basic::subs(ex): argument must be a list of equations"));
712 return subs(ls, lr, no_pattern);
715 /** Compare objects to establish canonical ordering.
716 * All compare functions return: -1 for *this less than other, 0 equal,
718 int basic::compare(const basic & other) const
720 unsigned hash_this = gethash();
721 unsigned hash_other = other.gethash();
723 if (hash_this<hash_other) return -1;
724 if (hash_this>hash_other) return 1;
726 unsigned typeid_this = tinfo();
727 unsigned typeid_other = other.tinfo();
729 if (typeid_this<typeid_other) {
730 // std::cout << "hash collision, different types: "
731 // << *this << " and " << other << std::endl;
732 // this->print(print_tree(std::cout));
733 // std::cout << " and ";
734 // other.print(print_tree(std::cout));
735 // std::cout << std::endl;
738 if (typeid_this>typeid_other) {
739 // std::cout << "hash collision, different types: "
740 // << *this << " and " << other << std::endl;
741 // this->print(print_tree(std::cout));
742 // std::cout << " and ";
743 // other.print(print_tree(std::cout));
744 // std::cout << std::endl;
748 GINAC_ASSERT(typeid(*this)==typeid(other));
750 // int cmpval = compare_same_type(other);
751 // if ((cmpval!=0) && (hash_this<0x80000000U)) {
752 // std::cout << "hash collision, same type: "
753 // << *this << " and " << other << std::endl;
754 // this->print(print_tree(std::cout));
755 // std::cout << " and ";
756 // other.print(print_tree(std::cout));
757 // std::cout << std::endl;
761 return compare_same_type(other);
764 /** Test for equality.
765 * This is only a quick test, meaning objects should be in the same domain.
766 * You might have to .expand(), .normal() objects first, depending on the
767 * domain of your computation, to get a more reliable answer.
769 * @see is_equal_same_type */
770 bool basic::is_equal(const basic & other) const
772 if (this->gethash()!=other.gethash())
774 if (this->tinfo()!=other.tinfo())
777 GINAC_ASSERT(typeid(*this)==typeid(other));
779 return is_equal_same_type(other);
784 /** Stop further evaluation.
786 * @see basic::eval */
787 const basic & basic::hold(void) const
789 return setflag(status_flags::evaluated);
792 /** Ensure the object may be modified without hurting others, throws if this
793 * is not the case. */
794 void basic::ensure_if_modifiable(void) const
796 if (this->refcount>1)
797 throw(std::runtime_error("cannot modify multiply referenced object"));
798 clearflag(status_flags::hash_calculated);
805 int max_recursion_level = 1024;