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
37 #include "relational.h"
38 #include "operators.h"
46 GINAC_IMPLEMENT_REGISTERED_CLASS_NO_CTORS(basic, void)
49 // default ctor, dtor, copy ctor, assignment operator and helpers
54 basic::basic(const basic & other) : tinfo_key(TINFO_basic), flags(0), refcount(0)
59 const basic & basic::operator=(const basic & other)
70 // none (all conditionally inlined)
76 // none (all conditionally inlined)
82 /** Construct object from archive_node. */
83 basic::basic(const archive_node &n, lst &sym_lst) : flags(0), refcount(0)
85 // Reconstruct tinfo_key from class name
86 std::string class_name;
87 if (n.find_string("class", class_name))
88 tinfo_key = find_tinfo_key(class_name);
90 throw (std::runtime_error("archive node contains no class name"));
93 /** Unarchive the object. */
94 DEFAULT_UNARCHIVE(basic)
96 /** Archive the object. */
97 void basic::archive(archive_node &n) const
99 n.add_string("class", class_name());
103 // new virtual functions which can be overridden by derived classes
108 /** Output to stream.
109 * @param c print context object that describes the output formatting
110 * @param level value that is used to identify the precedence or indentation
111 * level for placing parentheses and formatting */
112 void basic::print(const print_context & c, unsigned level) const
114 if (is_a<print_tree>(c)) {
116 c.s << std::string(level, ' ') << class_name()
117 << std::hex << ", hash=0x" << hashvalue << ", flags=0x" << flags << std::dec
118 << ", nops=" << nops()
120 for (size_t i=0; i<nops(); ++i)
121 op(i).print(c, level + static_cast<const print_tree &>(c).delta_indent);
124 c.s << "[" << class_name() << " object]";
127 /** Little wrapper around print to be called within a debugger.
128 * This is needed because you cannot call foo.print(cout) from within the
129 * debugger because it might not know what cout is. This method can be
130 * invoked with no argument and it will simply print to stdout.
133 * @see basic::dbgprinttree */
134 void basic::dbgprint(void) const
136 this->print(std::cerr);
137 std::cerr << std::endl;
140 /** Little wrapper around printtree to be called within a debugger.
142 * @see basic::dbgprint */
143 void basic::dbgprinttree(void) const
145 this->print(print_tree(std::cerr));
148 /** Return relative operator precedence (for parenthizing output). */
149 unsigned basic::precedence(void) const
154 /** Create a new copy of this on the heap. One can think of this as simulating
155 * a virtual copy constructor which is needed for instance by the refcounted
156 * construction of an ex from a basic. */
157 basic * basic::duplicate() const
159 return new basic(*this);
162 /** Information about the object.
164 * @see class info_flags */
165 bool basic::info(unsigned inf) const
167 // all possible properties are false for basic objects
171 /** Number of operands/members. */
172 size_t basic::nops() const
174 // iterating from 0 to nops() on atomic objects should be an empty loop,
175 // and accessing their elements is a range error. Container objects should
180 /** Return operand/member at position i. */
181 ex basic::op(size_t i) const
183 throw(std::range_error(std::string("basic::op(): ") + class_name() + std::string(" has no operands")));
186 /** Return modifyable operand/member at position i. */
187 ex & basic::let_op(size_t i)
189 ensure_if_modifiable();
190 throw(std::range_error(std::string("basic::let_op(): ") + class_name() + std::string(" has no operands")));
193 ex basic::operator[](const ex & index) const
195 if (is_exactly_a<numeric>(index))
196 return op(static_cast<size_t>(ex_to<numeric>(index).to_int()));
198 throw(std::invalid_argument(std::string("non-numeric indices not supported by ") + class_name()));
201 ex basic::operator[](size_t i) const
206 ex & basic::operator[](const ex & index)
208 if (is_exactly_a<numeric>(index))
209 return let_op(ex_to<numeric>(index).to_int());
211 throw(std::invalid_argument(std::string("non-numeric indices not supported by ") + class_name()));
214 ex & basic::operator[](size_t i)
219 /** Test for occurrence of a pattern. An object 'has' a pattern if it matches
220 * the pattern itself or one of the children 'has' it. As a consequence
221 * (according to the definition of children) given e=x+y+z, e.has(x) is true
222 * but e.has(x+y) is false. */
223 bool basic::has(const ex & pattern) const
226 if (match(pattern, repl_lst))
228 for (size_t i=0; i<nops(); i++)
229 if (op(i).has(pattern))
235 /** Construct new expression by applying the specified function to all
236 * sub-expressions (one level only, not recursively). */
237 ex basic::map(map_function & f) const
243 basic *copy = duplicate();
244 copy->setflag(status_flags::dynallocated);
245 copy->clearflag(status_flags::hash_calculated | status_flags::expanded);
247 for (size_t i=0; i<num; i++)
248 e.let_op(i) = f(e.op(i));
252 /** Return degree of highest power in object s. */
253 int basic::degree(const ex & s) const
255 return is_equal(ex_to<basic>(s)) ? 1 : 0;
258 /** Return degree of lowest power in object s. */
259 int basic::ldegree(const ex & s) const
261 return is_equal(ex_to<basic>(s)) ? 1 : 0;
264 /** Return coefficient of degree n in object s. */
265 ex basic::coeff(const ex & s, int n) const
267 if (is_equal(ex_to<basic>(s)))
268 return n==1 ? _ex1 : _ex0;
270 return n==0 ? *this : _ex0;
273 /** Sort expanded expression in terms of powers of some object(s).
274 * @param s object(s) to sort in
275 * @param distributed recursive or distributed form (only used when s is a list) */
276 ex basic::collect(const ex & s, bool distributed) const
281 // List of objects specified
285 return collect(s.op(0));
287 else if (distributed) {
289 // Get lower/upper degree of all symbols in list
290 size_t num = s.nops();
294 int cnt; // current degree, 'counter'
295 ex coeff; // coefficient for degree 'cnt'
297 sym_info *si = new sym_info[num];
299 for (size_t i=0; i<num; i++) {
301 si[i].ldeg = si[i].cnt = this->ldegree(si[i].sym);
302 si[i].deg = this->degree(si[i].sym);
303 c = si[i].coeff = c.coeff(si[i].sym, si[i].cnt);
308 // Calculate coeff*x1^c1*...*xn^cn
310 for (size_t i=0; i<num; i++) {
312 y *= power(si[i].sym, cnt);
314 x += y * si[num - 1].coeff;
316 // Increment counters
320 if (si[n].cnt <= si[n].deg) {
321 // Update coefficients
327 for (size_t i=n; i<num; i++)
328 c = si[i].coeff = c.coeff(si[i].sym, si[i].cnt);
333 si[n].cnt = si[n].ldeg;
344 size_t n = s.nops() - 1;
355 // Only one object specified
356 for (int n=this->ldegree(s); n<=this->degree(s); ++n)
357 x += this->coeff(s,n)*power(s,n);
360 // correct for lost fractional arguments and return
361 return x + (*this - x).expand();
364 /** Perform automatic non-interruptive term rewriting rules. */
365 ex basic::eval(int level) const
367 // There is nothing to do for basic objects:
371 /** Function object to be applied by basic::evalf(). */
372 struct evalf_map_function : public map_function {
374 evalf_map_function(int l) : level(l) {}
375 ex operator()(const ex & e) { return evalf(e, level); }
378 /** Evaluate object numerically. */
379 ex basic::evalf(int level) const
386 else if (level == -max_recursion_level)
387 throw(std::runtime_error("max recursion level reached"));
389 evalf_map_function map_evalf(level - 1);
390 return map(map_evalf);
395 /** Function object to be applied by basic::evalm(). */
396 struct evalm_map_function : public map_function {
397 ex operator()(const ex & e) { return evalm(e); }
400 /** Evaluate sums, products and integer powers of matrices. */
401 ex basic::evalm(void) const
406 return map(map_evalm);
409 /** Perform automatic symbolic evaluations on indexed expression that
410 * contains this object as the base expression. */
411 ex basic::eval_indexed(const basic & i) const
412 // this function can't take a "const ex & i" because that would result
413 // in an infinite eval() loop
415 // There is nothing to do for basic objects
419 /** Add two indexed expressions. They are guaranteed to be of class indexed
420 * (or a subclass) and their indices are compatible. This function is used
421 * internally by simplify_indexed().
423 * @param self First indexed expression; it's base object is *this
424 * @param other Second indexed expression
425 * @return sum of self and other
426 * @see ex::simplify_indexed() */
427 ex basic::add_indexed(const ex & self, const ex & other) const
432 /** Multiply an indexed expression with a scalar. This function is used
433 * internally by simplify_indexed().
435 * @param self Indexed expression; it's base object is *this
436 * @param other Numeric value
437 * @return product of self and other
438 * @see ex::simplify_indexed() */
439 ex basic::scalar_mul_indexed(const ex & self, const numeric & other) const
444 /** Try to contract two indexed expressions that appear in the same product.
445 * If a contraction exists, the function overwrites one or both of the
446 * expressions and returns true. Otherwise it returns false. It is
447 * guaranteed that both expressions are of class indexed (or a subclass)
448 * and that at least one dummy index has been found. This functions is
449 * used internally by simplify_indexed().
451 * @param self Pointer to first indexed expression; it's base object is *this
452 * @param other Pointer to second indexed expression
453 * @param v The complete vector of factors
454 * @return true if the contraction was successful, false otherwise
455 * @see ex::simplify_indexed() */
456 bool basic::contract_with(exvector::iterator self, exvector::iterator other, exvector & v) const
462 /** Check whether the expression matches a given pattern. For every wildcard
463 * object in the pattern, an expression of the form "wildcard == matching_expression"
464 * is added to repl_lst. */
465 bool basic::match(const ex & pattern, lst & repl_lst) const
468 Sweet sweet shapes, sweet sweet shapes,
469 That's the key thing, right right.
470 Feed feed face, feed feed shapes,
471 But who is the king tonight?
472 Who is the king tonight?
473 Pattern is the thing, the key thing-a-ling,
474 But who is the king of Pattern?
475 But who is the king, the king thing-a-ling,
476 Who is the king of Pattern?
477 Bog is the king, the king thing-a-ling,
478 Bog is the king of Pattern.
479 Ba bu-bu-bu-bu bu-bu-bu-bu-bu-bu bu-bu
480 Bog is the king of Pattern.
483 if (is_exactly_a<wildcard>(pattern)) {
485 // Wildcard matches anything, but check whether we already have found
486 // a match for that wildcard first (if so, it the earlier match must
487 // be the same expression)
488 for (lst::const_iterator it = repl_lst.begin(); it != repl_lst.end(); ++it) {
489 if (it->op(0).is_equal(pattern))
490 return is_equal(ex_to<basic>(it->op(1)));
492 repl_lst.append(pattern == *this);
497 // Expression must be of the same type as the pattern
498 if (tinfo() != ex_to<basic>(pattern).tinfo())
501 // Number of subexpressions must match
502 if (nops() != pattern.nops())
505 // No subexpressions? Then just compare the objects (there can't be
506 // wildcards in the pattern)
508 return is_equal_same_type(ex_to<basic>(pattern));
510 // Check whether attributes that are not subexpressions match
511 if (!match_same_type(ex_to<basic>(pattern)))
514 // Otherwise the subexpressions must match one-to-one
515 for (size_t i=0; i<nops(); i++)
516 if (!op(i).match(pattern.op(i), repl_lst))
519 // Looks similar enough, match found
524 /** Substitute a set of objects by arbitrary expressions. The ex returned
525 * will already be evaluated. */
526 ex basic::subs(const lst & ls, const lst & lr, unsigned options) const
528 GINAC_ASSERT(ls.nops() == lr.nops());
530 lst::const_iterator its, itr;
532 if (options & subs_options::no_pattern) {
533 for (its = ls.begin(), itr = lr.begin(); its != ls.end(); ++its, ++itr) {
534 if (is_equal(ex_to<basic>(*its)))
538 for (its = ls.begin(), itr = lr.begin(); its != ls.end(); ++its, ++itr) {
540 if (match(ex_to<basic>(*its), repl_lst))
541 return itr->subs(repl_lst, options | subs_options::no_pattern); // avoid infinite recursion when re-substituting the wildcards
548 /** Default interface of nth derivative ex::diff(s, n). It should be called
549 * instead of ::derivative(s) for first derivatives and for nth derivatives it
550 * just recurses down.
552 * @param s symbol to differentiate in
553 * @param nth order of differentiation
555 ex basic::diff(const symbol & s, unsigned nth) const
557 // trivial: zeroth derivative
561 // evaluate unevaluated *this before differentiating
562 if (!(flags & status_flags::evaluated))
563 return ex(*this).diff(s, nth);
565 ex ndiff = this->derivative(s);
566 while (!ndiff.is_zero() && // stop differentiating zeros
568 ndiff = ndiff.diff(s);
574 /** Return a vector containing the free indices of an expression. */
575 exvector basic::get_free_indices(void) const
577 return exvector(); // return an empty exvector
580 ex basic::eval_ncmul(const exvector & v) const
582 return hold_ncmul(v);
587 /** Function object to be applied by basic::derivative(). */
588 struct derivative_map_function : public map_function {
590 derivative_map_function(const symbol &sym) : s(sym) {}
591 ex operator()(const ex & e) { return diff(e, s); }
594 /** Default implementation of ex::diff(). It maps the operation on the
595 * operands (or returns 0 when the object has no operands).
598 ex basic::derivative(const symbol & s) const
603 derivative_map_function map_derivative(s);
604 return map(map_derivative);
608 /** Returns order relation between two objects of same type. This needs to be
609 * implemented by each class. It may never return anything else than 0,
610 * signalling equality, or +1 and -1 signalling inequality and determining
611 * the canonical ordering. (Perl hackers will wonder why C++ doesn't feature
612 * the spaceship operator <=> for denoting just this.) */
613 int basic::compare_same_type(const basic & other) const
615 return compare_pointers(this, &other);
618 /** Returns true if two objects of same type are equal. Normally needs
619 * not be reimplemented as long as it wasn't overwritten by some parent
620 * class, since it just calls compare_same_type(). The reason why this
621 * function exists is that sometimes it is easier to determine equality
622 * than an order relation and then it can be overridden. */
623 bool basic::is_equal_same_type(const basic & other) const
625 return compare_same_type(other)==0;
628 /** Returns true if the attributes of two objects are similar enough for
629 * a match. This function must not match subexpressions (this is already
630 * done by basic::match()). Only attributes not accessible by op() should
631 * be compared. This is also the reason why this function doesn't take the
632 * wildcard replacement list from match() as an argument: only subexpressions
633 * are subject to wildcard matches. Also, this function only needs to be
634 * implemented for container classes because is_equal_same_type() is
635 * automatically used instead of match_same_type() if nops() == 0.
637 * @see basic::match */
638 bool basic::match_same_type(const basic & other) const
640 // The default is to only consider subexpressions, but not any other
645 unsigned basic::return_type(void) const
647 return return_types::commutative;
650 unsigned basic::return_type_tinfo(void) const
655 /** Compute the hash value of an object and if it makes sense to store it in
656 * the objects status_flags, do so. The method inherited from class basic
657 * computes a hash value based on the type and hash values of possible
658 * members. For this reason it is well suited for container classes but
659 * atomic classes should override this implementation because otherwise they
660 * would all end up with the same hashvalue. */
661 unsigned basic::calchash(void) const
663 unsigned v = golden_ratio_hash(tinfo());
664 for (size_t i=0; i<nops(); i++) {
666 v ^= (const_cast<basic *>(this))->op(i).gethash();
669 // store calculated hash value only if object is already evaluated
670 if (flags & status_flags::evaluated) {
671 setflag(status_flags::hash_calculated);
678 /** Function object to be applied by basic::expand(). */
679 struct expand_map_function : public map_function {
681 expand_map_function(unsigned o) : options(o) {}
682 ex operator()(const ex & e) { return expand(e, options); }
685 /** Expand expression, i.e. multiply it out and return the result as a new
687 ex basic::expand(unsigned options) const
690 return (options == 0) ? setflag(status_flags::expanded) : *this;
692 expand_map_function map_expand(options);
693 return ex_to<basic>(map(map_expand)).setflag(options == 0 ? status_flags::expanded : 0);
699 // non-virtual functions in this class
704 /** Substitute objects in an expression (syntactic substitution) and return
705 * the result as a new expression. There are two valid types of
706 * replacement arguments: 1) a relational like object==ex and 2) a list of
707 * relationals lst(object1==ex1,object2==ex2,...), which is converted to
708 * subs(lst(object1,object2,...),lst(ex1,ex2,...)). */
709 ex basic::subs(const ex & e, unsigned options) const
711 if (e.info(info_flags::relation_equal)) {
712 return subs(lst(e.lhs()), lst(e.rhs()), options);
714 if (!e.info(info_flags::list)) {
715 throw(std::invalid_argument("basic::subs(ex): argument must be a list"));
718 // Split list into two
721 GINAC_ASSERT(is_a<lst>(e));
722 for (lst::const_iterator it = ex_to<lst>(e).begin(); it != ex_to<lst>(e).end(); ++it) {
724 if (!r.info(info_flags::relation_equal)) {
725 throw(std::invalid_argument("basic::subs(ex): argument must be a list of equations"));
727 const ex & s = r.op(0);
731 // Search for products and powers in the expressions to be substituted
732 // (for an optimization in expairseq::subs())
733 if (is_exactly_a<mul>(s) || is_exactly_a<power>(s))
734 options |= subs_options::pattern_is_product;
736 if (!(options & subs_options::pattern_is_product))
737 options |= subs_options::pattern_is_not_product;
739 return subs(ls, lr, options);
742 /** Compare objects syntactically to establish canonical ordering.
743 * All compare functions return: -1 for *this less than other, 0 equal,
745 int basic::compare(const basic & other) const
747 const unsigned hash_this = gethash();
748 const unsigned hash_other = other.gethash();
749 if (hash_this<hash_other) return -1;
750 if (hash_this>hash_other) return 1;
752 const unsigned typeid_this = tinfo();
753 const unsigned typeid_other = other.tinfo();
754 if (typeid_this==typeid_other) {
755 GINAC_ASSERT(typeid(*this)==typeid(other));
756 // int cmpval = compare_same_type(other);
758 // std::cout << "hash collision, same type: "
759 // << *this << " and " << other << std::endl;
760 // this->print(print_tree(std::cout));
761 // std::cout << " and ";
762 // other.print(print_tree(std::cout));
763 // std::cout << std::endl;
766 return compare_same_type(other);
768 // std::cout << "hash collision, different types: "
769 // << *this << " and " << other << std::endl;
770 // this->print(print_tree(std::cout));
771 // std::cout << " and ";
772 // other.print(print_tree(std::cout));
773 // std::cout << std::endl;
774 return (typeid_this<typeid_other ? -1 : 1);
778 /** Test for syntactic equality.
779 * This is only a quick test, meaning objects should be in the same domain.
780 * You might have to .expand(), .normal() objects first, depending on the
781 * domain of your computation, to get a more reliable answer.
783 * @see is_equal_same_type */
784 bool basic::is_equal(const basic & other) const
786 if (this->gethash()!=other.gethash())
788 if (this->tinfo()!=other.tinfo())
791 GINAC_ASSERT(typeid(*this)==typeid(other));
793 return is_equal_same_type(other);
798 /** Stop further evaluation.
800 * @see basic::eval */
801 const basic & basic::hold(void) const
803 return setflag(status_flags::evaluated);
806 /** Ensure the object may be modified without hurting others, throws if this
807 * is not the case. */
808 void basic::ensure_if_modifiable(void) const
811 throw(std::runtime_error("cannot modify multiply referenced object"));
812 clearflag(status_flags::hash_calculated);
819 int max_recursion_level = 1024;