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)
54 debugmsg("basic copy ctor", LOGLEVEL_CONSTRUCT);
58 const basic & basic::operator=(const basic & other)
60 debugmsg("basic operator=", LOGLEVEL_ASSIGNMENT);
70 // none (all conditionally inlined)
76 // none (all conditionally inlined)
82 /** Construct object from archive_node. */
83 basic::basic(const archive_node &n, const lst &sym_lst) : flags(0), refcount(0)
85 debugmsg("basic ctor from archive_node", LOGLEVEL_CONSTRUCT);
87 // Reconstruct tinfo_key from class name
88 std::string class_name;
89 if (n.find_string("class", class_name))
90 tinfo_key = find_tinfo_key(class_name);
92 throw (std::runtime_error("archive node contains no class name"));
95 /** Unarchive the object. */
96 DEFAULT_UNARCHIVE(basic)
98 /** Archive the object. */
99 void basic::archive(archive_node &n) const
101 n.add_string("class", class_name());
105 // functions overriding virtual functions from bases classes
111 // new virtual functions which can be overridden by derived classes
116 /** Output to stream.
117 * @param c print context object that describes the output formatting
118 * @param level value that is used to identify the precedence or indentation
119 * level for placing parentheses and formatting */
120 void basic::print(const print_context & c, unsigned level) const
122 debugmsg("basic print", LOGLEVEL_PRINT);
124 if (is_of_type(c, print_tree)) {
126 c.s << std::string(level, ' ') << class_name()
127 << std::hex << ", hash=0x" << hashvalue << ", flags=0x" << flags << std::dec
128 << ", nops=" << nops()
130 for (unsigned i=0; i<nops(); ++i)
131 op(i).print(c, level + static_cast<const print_tree &>(c).delta_indent);
134 c.s << "[" << class_name() << " object]";
137 /** Little wrapper around print to be called within a debugger.
138 * This is needed because you cannot call foo.print(cout) from within the
139 * debugger because it might not know what cout is. This method can be
140 * invoked with no argument and it will simply print to stdout.
142 * @see basic::print */
143 void basic::dbgprint(void) const
145 this->print(std::cerr);
146 std::cerr << std::endl;
149 /** Little wrapper around printtree to be called within a debugger.
151 * @see basic::dbgprint
152 * @see basic::printtree */
153 void basic::dbgprinttree(void) const
155 this->print(print_tree(std::cerr));
158 /** Return relative operator precedence (for parenthizing output). */
159 unsigned basic::precedence(void) const
164 /** Create a new copy of this on the heap. One can think of this as simulating
165 * a virtual copy constructor which is needed for instance by the refcounted
166 * construction of an ex from a basic. */
167 basic * basic::duplicate() const
169 debugmsg("basic duplicate",LOGLEVEL_DUPLICATE);
170 return new basic(*this);
173 /** Information about the object.
175 * @see class info_flags */
176 bool basic::info(unsigned inf) const
178 // all possible properties are false for basic objects
182 /** Number of operands/members. */
183 unsigned basic::nops() const
185 // iterating from 0 to nops() on atomic objects should be an empty loop,
186 // and accessing their elements is a range error. Container objects should
191 /** Return operand/member at position i. */
192 ex basic::op(int i) const
194 return (const_cast<basic *>(this))->let_op(i);
197 /** Return modifyable operand/member at position i. */
198 ex & basic::let_op(int i)
200 throw(std::out_of_range("op() out of range"));
203 ex basic::operator[](const ex & index) const
205 if (is_exactly_of_type(*index.bp,numeric))
206 return op(static_cast<const numeric &>(*index.bp).to_int());
208 throw(std::invalid_argument("non-numeric indices not supported by this type"));
211 ex basic::operator[](int i) const
216 /** Test for occurrence of a pattern. An object 'has' a pattern if it matches
217 * the pattern itself or one of the children 'has' it. As a consequence
218 * (according to the definition of children) given e=x+y+z, e.has(x) is true
219 * but e.has(x+y) is false. */
220 bool basic::has(const ex & pattern) const
222 GINAC_ASSERT(pattern.bp!=0);
224 if (match(*pattern.bp, repl_lst))
226 for (unsigned i=0; i<nops(); i++)
227 if (op(i).has(pattern))
233 /** Construct new expression by applying the specified function to all
234 * sub-expressions (one level only, not recursively). */
235 ex basic::map(map_function & f) const
237 unsigned num = nops();
241 basic *copy = duplicate();
242 copy->setflag(status_flags::dynallocated);
243 copy->clearflag(status_flags::hash_calculated | status_flags::expanded);
245 for (unsigned i=0; i<num; i++)
246 e.let_op(i) = f(e.op(i));
250 /** Return degree of highest power in object s. */
251 int basic::degree(const ex & s) const
256 /** Return degree of lowest power in object s. */
257 int basic::ldegree(const ex & s) const
262 /** Return coefficient of degree n in object s. */
263 ex basic::coeff(const ex & s, int n) const
265 return n==0 ? *this : _ex0();
268 /** Sort expanded expression in terms of powers of some object(s).
269 * @param s object(s) to sort in
270 * @param distributed recursive or distributed form (only used when s is a list) */
271 ex basic::collect(const ex & s, bool distributed) const
274 if (is_ex_of_type(s, lst)) {
276 // List of objects specified
280 return collect(s.op(0));
282 else if (distributed) {
284 // Get lower/upper degree of all symbols in list
289 int cnt; // current degree, 'counter'
290 ex coeff; // coefficient for degree 'cnt'
292 sym_info *si = new sym_info[num];
294 for (int i=0; i<num; i++) {
296 si[i].ldeg = si[i].cnt = this->ldegree(si[i].sym);
297 si[i].deg = this->degree(si[i].sym);
298 c = si[i].coeff = c.coeff(si[i].sym, si[i].cnt);
303 // Calculate coeff*x1^c1*...*xn^cn
305 for (int i=0; i<num; i++) {
307 y *= power(si[i].sym, cnt);
309 x += y * si[num - 1].coeff;
311 // Increment counters
315 if (si[n].cnt <= si[n].deg) {
316 // Update coefficients
322 for (int i=n; i<num; i++)
323 c = si[i].coeff = c.coeff(si[i].sym, si[i].cnt);
328 si[n].cnt = si[n].ldeg;
339 for (int n=s.nops()-1; n>=0; n--)
345 // Only one object specified
346 for (int n=this->ldegree(s); n<=this->degree(s); ++n)
347 x += this->coeff(s,n)*power(s,n);
350 // correct for lost fractional arguments and return
351 return x + (*this - x).expand();
354 /** Perform automatic non-interruptive symbolic evaluation on expression. */
355 ex basic::eval(int level) const
357 // There is nothing to do for basic objects:
361 /** Function object to be applied by basic::evalf(). */
362 struct evalf_map_function : public map_function {
364 evalf_map_function(int l) : level(l) {}
365 ex operator()(const ex & e) { return evalf(e, level); }
368 /** Evaluate object numerically. */
369 ex basic::evalf(int level) const
376 else if (level == -max_recursion_level)
377 throw(std::runtime_error("max recursion level reached"));
379 evalf_map_function map_evalf(level - 1);
380 return map(map_evalf);
385 /** Function object to be applied by basic::evalm(). */
386 struct evalm_map_function : public map_function {
387 ex operator()(const ex & e) { return evalm(e); }
390 /** Evaluate sums, products and integer powers of matrices. */
391 ex basic::evalm(void) const
396 return map(map_evalm);
399 /** Perform automatic symbolic evaluations on indexed expression that
400 * contains this object as the base expression. */
401 ex basic::eval_indexed(const basic & i) const
402 // this function can't take a "const ex & i" because that would result
403 // in an infinite eval() loop
405 // There is nothing to do for basic objects
409 /** Add two indexed expressions. They are guaranteed to be of class indexed
410 * (or a subclass) and their indices are compatible. This function is used
411 * internally by simplify_indexed().
413 * @param self First indexed expression; it's base object is *this
414 * @param other Second indexed expression
415 * @return sum of self and other
416 * @see ex::simplify_indexed() */
417 ex basic::add_indexed(const ex & self, const ex & other) const
422 /** Multiply an indexed expression with a scalar. This function is used
423 * internally by simplify_indexed().
425 * @param self Indexed expression; it's base object is *this
426 * @param other Numeric value
427 * @return product of self and other
428 * @see ex::simplify_indexed() */
429 ex basic::scalar_mul_indexed(const ex & self, const numeric & other) const
434 /** Try to contract two indexed expressions that appear in the same product.
435 * If a contraction exists, the function overwrites one or both of the
436 * expressions and returns true. Otherwise it returns false. It is
437 * guaranteed that both expressions are of class indexed (or a subclass)
438 * and that at least one dummy index has been found. This functions is
439 * used internally by simplify_indexed().
441 * @param self Pointer to first indexed expression; it's base object is *this
442 * @param other Pointer to second indexed expression
443 * @param v The complete vector of factors
444 * @return true if the contraction was successful, false otherwise
445 * @see ex::simplify_indexed() */
446 bool basic::contract_with(exvector::iterator self, exvector::iterator other, exvector & v) const
452 /** Check whether the expression matches a given pattern. For every wildcard
453 * object in the pattern, an expression of the form "wildcard == matching_expression"
454 * is added to repl_lst. */
455 bool basic::match(const ex & pattern, lst & repl_lst) const
458 Sweet sweet shapes, sweet sweet shapes,
459 That's the key thing, right right.
460 Feed feed face, feed feed shapes,
461 But who is the king tonight?
462 Who is the king tonight?
463 Pattern is the thing, the key thing-a-ling,
464 But who is the king of Pattern?
465 But who is the king, the king thing-a-ling,
466 Who is the king of Pattern?
467 Bog is the king, the king thing-a-ling,
468 Bog is the king of Pattern.
469 Ba bu-bu-bu-bu bu-bu-bu-bu-bu-bu bu-bu
470 Bog is the king of Pattern.
473 if (is_ex_exactly_of_type(pattern, wildcard)) {
475 // Wildcard matches anything, but check whether we already have found
476 // a match for that wildcard first (if so, it the earlier match must
477 // be the same expression)
478 for (unsigned i=0; i<repl_lst.nops(); i++) {
479 if (repl_lst.op(i).op(0).is_equal(pattern))
480 return is_equal(*repl_lst.op(i).op(1).bp);
482 repl_lst.append(pattern == *this);
487 // Expression must be of the same type as the pattern
488 if (tinfo() != pattern.bp->tinfo())
491 // Number of subexpressions must match
492 if (nops() != pattern.nops())
495 // No subexpressions? Then just compare the objects (there can't be
496 // wildcards in the pattern)
498 return is_equal_same_type(*pattern.bp);
500 // Check whether attributes that are not subexpressions match
501 if (!match_same_type(*pattern.bp))
504 // Otherwise the subexpressions must match one-to-one
505 for (unsigned i=0; i<nops(); i++)
506 if (!op(i).match(pattern.op(i), repl_lst))
509 // Looks similar enough, match found
514 /** Substitute a set of objects by arbitrary expressions. The ex returned
515 * will already be evaluated. */
516 ex basic::subs(const lst & ls, const lst & lr, bool no_pattern) const
518 GINAC_ASSERT(ls.nops() == lr.nops());
521 for (unsigned i=0; i<ls.nops(); i++) {
522 if (is_equal(*ls.op(i).bp))
526 for (unsigned i=0; i<ls.nops(); i++) {
528 if (match(*ls.op(i).bp, repl_lst))
529 return lr.op(i).bp->subs(repl_lst, true); // avoid infinite recursion when re-substituting the wildcards
536 /** Default interface of nth derivative ex::diff(s, n). It should be called
537 * instead of ::derivative(s) for first derivatives and for nth derivatives it
538 * just recurses down.
540 * @param s symbol to differentiate in
541 * @param nth order of differentiation
543 ex basic::diff(const symbol & s, unsigned nth) const
545 // trivial: zeroth derivative
549 // evaluate unevaluated *this before differentiating
550 if (!(flags & status_flags::evaluated))
551 return ex(*this).diff(s, nth);
553 ex ndiff = this->derivative(s);
554 while (!ndiff.is_zero() && // stop differentiating zeros
556 ndiff = ndiff.diff(s);
562 /** Return a vector containing the free indices of an expression. */
563 exvector basic::get_free_indices(void) const
565 return exvector(); // return an empty exvector
568 ex basic::simplify_ncmul(const exvector & v) const
570 return simplified_ncmul(v);
575 /** Function object to be applied by basic::derivative(). */
576 struct derivative_map_function : public map_function {
578 derivative_map_function(const symbol &sym) : s(sym) {}
579 ex operator()(const ex & e) { return diff(e, s); }
582 /** Default implementation of ex::diff(). It maps the operation on the
583 * operands (or returns 0 when the object has no operands).
586 ex basic::derivative(const symbol & s) const
591 derivative_map_function map_derivative(s);
592 return map(map_derivative);
596 /** Returns order relation between two objects of same type. This needs to be
597 * implemented by each class. It may never return anything else than 0,
598 * signalling equality, or +1 and -1 signalling inequality and determining
599 * the canonical ordering. (Perl hackers will wonder why C++ doesn't feature
600 * the spaceship operator <=> for denoting just this.) */
601 int basic::compare_same_type(const basic & other) const
603 return compare_pointers(this, &other);
606 /** Returns true if two objects of same type are equal. Normally needs
607 * not be reimplemented as long as it wasn't overwritten by some parent
608 * class, since it just calls compare_same_type(). The reason why this
609 * function exists is that sometimes it is easier to determine equality
610 * than an order relation and then it can be overridden. */
611 bool basic::is_equal_same_type(const basic & other) const
613 return compare_same_type(other)==0;
616 /** Returns true if the attributes of two objects are similar enough for
617 * a match. This function must not match subexpressions (this is already
618 * done by basic::match()). Only attributes not accessible by op() should
619 * be compared. This is also the reason why this function doesn't take the
620 * wildcard replacement list from match() as an argument: only subexpressions
621 * are subject to wildcard matches. Also, this function only needs to be
622 * implemented for container classes because is_equal_same_type() is
623 * automatically used instead of match_same_type() if nops() == 0.
625 * @see basic::match */
626 bool basic::match_same_type(const basic & other) const
628 // The default is to only consider subexpressions, but not any other
633 unsigned basic::return_type(void) const
635 return return_types::commutative;
638 unsigned basic::return_type_tinfo(void) const
643 /** Compute the hash value of an object and if it makes sense to store it in
644 * the objects status_flags, do so. The method inherited from class basic
645 * computes a hash value based on the type and hash values of possible
646 * members. For this reason it is well suited for container classes but
647 * atomic classes should override this implementation because otherwise they
648 * would all end up with the same hashvalue. */
649 unsigned basic::calchash(void) const
651 unsigned v = golden_ratio_hash(tinfo());
652 for (unsigned i=0; i<nops(); i++) {
653 v = rotate_left_31(v);
654 v ^= (const_cast<basic *>(this))->op(i).gethash();
657 // mask out numeric hashes:
660 // store calculated hash value only if object is already evaluated
661 if (flags & status_flags::evaluated) {
662 setflag(status_flags::hash_calculated);
669 /** Function object to be applied by basic::expand(). */
670 struct expand_map_function : public map_function {
672 expand_map_function(unsigned o) : options(o) {}
673 ex operator()(const ex & e) { return expand(e, options); }
676 /** Expand expression, i.e. multiply it out and return the result as a new
678 ex basic::expand(unsigned options) const
681 return (options == 0) ? setflag(status_flags::expanded) : *this;
683 expand_map_function map_expand(options);
684 return map(map_expand).bp->setflag(options == 0 ? status_flags::expanded : 0);
690 // non-virtual functions in this class
695 /** Substitute objects in an expression (syntactic substitution) and return
696 * the result as a new expression. There are two valid types of
697 * replacement arguments: 1) a relational like object==ex and 2) a list of
698 * relationals lst(object1==ex1,object2==ex2,...), which is converted to
699 * subs(lst(object1,object2,...),lst(ex1,ex2,...)). */
700 ex basic::subs(const ex & e, bool no_pattern) const
702 if (e.info(info_flags::relation_equal)) {
703 return subs(lst(e), no_pattern);
705 if (!e.info(info_flags::list)) {
706 throw(std::invalid_argument("basic::subs(ex): argument must be a list"));
710 for (unsigned i=0; i<e.nops(); i++) {
712 if (!r.info(info_flags::relation_equal)) {
713 throw(std::invalid_argument("basic::subs(ex): argument must be a list of equations"));
718 return subs(ls, lr, no_pattern);
721 /** Compare objects to establish canonical ordering.
722 * All compare functions return: -1 for *this less than other, 0 equal,
724 int basic::compare(const basic & other) const
726 unsigned hash_this = gethash();
727 unsigned hash_other = other.gethash();
729 if (hash_this<hash_other) return -1;
730 if (hash_this>hash_other) return 1;
732 unsigned typeid_this = tinfo();
733 unsigned typeid_other = other.tinfo();
735 if (typeid_this<typeid_other) {
736 // std::cout << "hash collision, different types: "
737 // << *this << " and " << other << std::endl;
738 // this->print(print_tree(std::cout));
739 // std::cout << " and ";
740 // other.print(print_tree(std::cout));
741 // std::cout << std::endl;
744 if (typeid_this>typeid_other) {
745 // std::cout << "hash collision, different types: "
746 // << *this << " and " << other << std::endl;
747 // this->print(print_tree(std::cout));
748 // std::cout << " and ";
749 // other.print(print_tree(std::cout));
750 // std::cout << std::endl;
754 GINAC_ASSERT(typeid(*this)==typeid(other));
756 // int cmpval = compare_same_type(other);
757 // if ((cmpval!=0) && (hash_this<0x80000000U)) {
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;
767 return compare_same_type(other);
770 /** Test for equality.
771 * This is only a quick test, meaning objects should be in the same domain.
772 * You might have to .expand(), .normal() objects first, depending on the
773 * domain of your computation, to get a more reliable answer.
775 * @see is_equal_same_type */
776 bool basic::is_equal(const basic & other) const
778 if (this->gethash()!=other.gethash())
780 if (this->tinfo()!=other.tinfo())
783 GINAC_ASSERT(typeid(*this)==typeid(other));
785 return is_equal_same_type(other);
790 /** Stop further evaluation.
792 * @see basic::eval */
793 const basic & basic::hold(void) const
795 return setflag(status_flags::evaluated);
798 /** Ensure the object may be modified without hurting others, throws if this
799 * is not the case. */
800 void basic::ensure_if_modifiable(void) const
802 if (this->refcount>1)
803 throw(std::runtime_error("cannot modify multiply referenced object"));
804 clearflag(status_flags::hash_calculated);
811 int max_recursion_level = 1024;