ginac/basic.cpp

   1 /** @file basic.cpp
   2  *
   3  *  Implementation of GiNaC's ABC. */
   4
   5 /*
   6  *  GiNaC Copyright (C) 1999-2001 Johannes Gutenberg University Mainz, Germany
   7  *
   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.
  12  *
  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.
  17  *
  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
  21  */
  22
  23 #include <iostream>
  24 #include <stdexcept>
  25 #ifdef DO_GINAC_ASSERT
  26 #  include <typeinfo>
  27 #endif
  28
  29 #include "basic.h"
  30 #include "ex.h"
  31 #include "numeric.h"
  32 #include "power.h"
  33 #include "symbol.h"
  34 #include "lst.h"
  35 #include "ncmul.h"
  36 #include "relational.h"
  37 #include "print.h"
  38 #include "archive.h"
  39 #include "utils.h"
  40 #include "debugmsg.h"
  41
  42 namespace GiNaC {
  43
  44 GINAC_IMPLEMENT_REGISTERED_CLASS_NO_CTORS(basic, void)
  45
  46 //////////
  47 // default ctor, dtor, copy ctor assignment operator and helpers
  48 //////////
  49
  50 // public
  51
  52 basic::basic(const basic & other) : tinfo_key(TINFO_basic), flags(0), refcount(0)
  53 {
  54         debugmsg("basic copy ctor", LOGLEVEL_CONSTRUCT);
  55         copy(other);
  56 }
  57
  58 const basic & basic::operator=(const basic & other)
  59 {
  60         debugmsg("basic operator=", LOGLEVEL_ASSIGNMENT);
  61         if (this != &other) {
  62                 destroy(true);
  63                 copy(other);
  64         }
  65         return *this;
  66 }
  67
  68 // protected
  69
  70 // none (all conditionally inlined)
  71
  72 //////////
  73 // other ctors
  74 //////////
  75
  76 // none (all conditionally inlined)
  77
  78 //////////
  79 // archiving
  80 //////////
  81
  82 /** Construct object from archive_node. */
  83 basic::basic(const archive_node &n, const lst &sym_lst) : flags(0), refcount(0)
  84 {
  85         debugmsg("basic ctor from archive_node", LOGLEVEL_CONSTRUCT);
  86
  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);
  91         else
  92                 throw (std::runtime_error("archive node contains no class name"));
  93 }
  94
  95 /** Unarchive the object. */
  96 DEFAULT_UNARCHIVE(basic)
  97
  98 /** Archive the object. */
  99 void basic::archive(archive_node &n) const
 100 {
 101         n.add_string("class", class_name());
 102 }
 103
 104 //////////
 105 // functions overriding virtual functions from bases classes
 106 //////////
 107
 108 // none
 109
 110 //////////
 111 // new virtual functions which can be overridden by derived classes
 112 //////////
 113
 114 // public
 115
 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
 121 {
 122         debugmsg("basic print", LOGLEVEL_PRINT);
 123
 124         if (is_of_type(c, print_tree)) {
 125
 126                 c.s << std::string(level, ' ') << class_name()
 127                     << std::hex << ", hash=0x" << hashvalue << ", flags=0x" << flags << std::dec
 128                     << ", nops=" << nops()
 129                     << std::endl;
 130                 for (unsigned i=0; i<nops(); ++i)
 131                         op(i).print(c, level + static_cast<const print_tree &>(c).delta_indent);
 132
 133         } else
 134                 c.s << "[" << class_name() << " object]";
 135 }
 136
 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.
 141  *
 142  *  @see basic::print */
 143 void basic::dbgprint(void) const
 144 {
 145         this->print(std::cerr);
 146         std::cerr << std::endl;
 147 }
 148
 149 /** Little wrapper around printtree to be called within a debugger.
 150  *
 151  *  @see basic::dbgprint
 152  *  @see basic::printtree */
 153 void basic::dbgprinttree(void) const
 154 {
 155         this->print(print_tree(std::cerr));
 156 }
 157
 158 /** Return relative operator precedence (for parenthizing output). */
 159 unsigned basic::precedence(void) const
 160 {
 161         return 70;
 162 }
 163
 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
 168 {
 169         debugmsg("basic duplicate",LOGLEVEL_DUPLICATE);
 170         return new basic(*this);
 171 }
 172
 173 /** Information about the object.
 174  *
 175  *  @see class info_flags */
 176 bool basic::info(unsigned inf) const
 177 {
 178         // all possible properties are false for basic objects
 179         return false;
 180 }
 181
 182 /** Number of operands/members. */
 183 unsigned basic::nops() const
 184 {
 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
 187         // override this.
 188         return 0;
 189 }
 190
 191 /** Return operand/member at position i. */
 192 ex basic::op(int i) const
 193 {
 194         return (const_cast<basic *>(this))->let_op(i);
 195 }
 196
 197 /** Return modifyable operand/member at position i. */
 198 ex & basic::let_op(int i)
 199 {
 200         throw(std::out_of_range("op() out of range"));
 201 }
 202
 203 ex basic::operator[](const ex & index) const
 204 {
 205         if (is_exactly_of_type(*index.bp,numeric))
 206                 return op(static_cast<const numeric &>(*index.bp).to_int());
 207
 208         throw(std::invalid_argument("non-numeric indices not supported by this type"));
 209 }
 210
 211 ex basic::operator[](int i) const
 212 {
 213         return op(i);
 214 }
 215
 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
 221 {
 222         GINAC_ASSERT(pattern.bp!=0);
 223         lst repl_lst;
 224         if (match(*pattern.bp, repl_lst))
 225                 return true;
 226         for (unsigned i=0; i<nops(); i++)
 227                 if (op(i).has(pattern))
 228                         return true;
 229
 230         return false;
 231 }
 232
 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
 236 {
 237         unsigned num = nops();
 238         if (num == 0)
 239                 return *this;
 240
 241         basic *copy = duplicate();
 242         copy->setflag(status_flags::dynallocated);
 243         copy->clearflag(status_flags::hash_calculated | status_flags::expanded);
 244         ex e(*copy);
 245         for (unsigned i=0; i<num; i++)
 246                 e.let_op(i) = f(e.op(i));
 247         return e.eval();
 248 }
 249
 250 /** Return degree of highest power in object s. */
 251 int basic::degree(const ex & s) const
 252 {
 253         return 0;
 254 }
 255
 256 /** Return degree of lowest power in object s. */
 257 int basic::ldegree(const ex & s) const
 258 {
 259         return 0;
 260 }
 261
 262 /** Return coefficient of degree n in object s. */
 263 ex basic::coeff(const ex & s, int n) const
 264 {
 265         return n==0 ? *this : _ex0();
 266 }
 267
 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
 272 {
 273         ex x;
 274         if (is_ex_of_type(s, lst)) {
 275
 276                 // List of objects specified
 277                 if (s.nops() == 0)
 278                         return *this;
 279                 if (s.nops() == 1)
 280                         return collect(s.op(0));
 281
 282                 else if (distributed) {
 283
 284                         // Get lower/upper degree of all symbols in list
 285                         int num = s.nops();
 286                         struct sym_info {
 287                                 ex sym;
 288                                 int ldeg, deg;
 289                                 int cnt;  // current degree, 'counter'
 290                                 ex coeff; // coefficient for degree 'cnt'
 291                         };
 292                         sym_info *si = new sym_info[num];
 293                         ex c = *this;
 294                         for (int i=0; i<num; i++) {
 295                                 si[i].sym = s.op(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);
 299                         }
 300
 301                         while (true) {
 302
 303                                 // Calculate coeff*x1^c1*...*xn^cn
 304                                 ex y = _ex1();
 305                                 for (int i=0; i<num; i++) {
 306                                         int cnt = si[i].cnt;
 307                                         y *= power(si[i].sym, cnt);
 308                                 }
 309                                 x += y * si[num - 1].coeff;
 310
 311                                 // Increment counters
 312                                 int n = num - 1;
 313                                 while (true) {
 314                                         ++si[n].cnt;
 315                                         if (si[n].cnt <= si[n].deg) {
 316                                                 // Update coefficients
 317                                                 ex c;
 318                                                 if (n == 0)
 319                                                         c = *this;
 320                                                 else
 321                                                         c = si[n - 1].coeff;
 322                                                 for (int i=n; i<num; i++)
 323                                                         c = si[i].coeff = c.coeff(si[i].sym, si[i].cnt);
 324                                                 break;
 325                                         }
 326                                         if (n == 0)
 327                                                 goto done;
 328                                         si[n].cnt = si[n].ldeg;
 329                                         n--;
 330                                 }
 331                         }
 332
 333 done:           delete[] si;
 334
 335                 } else {
 336
 337                         // Recursive form
 338                         x = *this;
 339                         for (int n=s.nops()-1; n>=0; n--)
 340                                 x = x.collect(s[n]);
 341                 }
 342
 343         } else {
 344
 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);
 348         }
 349
 350         // correct for lost fractional arguments and return
 351         return x + (*this - x).expand();
 352 }
 353
 354 /** Perform automatic non-interruptive symbolic evaluation on expression. */
 355 ex basic::eval(int level) const
 356 {
 357         // There is nothing to do for basic objects:
 358         return this->hold();
 359 }
 360
 361 /** Function object to be applied by basic::evalf(). */
 362 struct evalf_map_function : public map_function {
 363         int level;
 364         evalf_map_function(int l) : level(l) {}
 365         ex operator()(const ex & e) { return evalf(e, level); }
 366 };
 367
 368 /** Evaluate object numerically. */
 369 ex basic::evalf(int level) const
 370 {
 371         if (nops() == 0)
 372                 return *this;
 373         else {
 374                 if (level == 1)
 375                         return *this;
 376                 else if (level == -max_recursion_level)
 377                         throw(std::runtime_error("max recursion level reached"));
 378                 else {
 379                         evalf_map_function map_evalf(level - 1);
 380                         return map(map_evalf);
 381                 }
 382         }
 383 }
 384
 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); }
 388 } map_evalm;
 389
 390 /** Evaluate sums, products and integer powers of matrices. */
 391 ex basic::evalm(void) const
 392 {
 393         if (nops() == 0)
 394                 return *this;
 395         else
 396                 return map(map_evalm);
 397 }
 398
 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
 404 {
 405         // There is nothing to do for basic objects
 406         return i.hold();
 407 }
 408
 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().
 412  *
 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
 418 {
 419         return self + other;
 420 }
 421
 422 /** Multiply an indexed expression with a scalar. This function is used
 423  *  internally by simplify_indexed().
 424  *
 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
 430 {
 431         return self * other;
 432 }
 433
 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().
 440  *
 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
 447 {
 448         // Do nothing
 449         return false;
 450 }
 451
 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
 456 {
 457 /*
 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.
 471 */
 472
 473         if (is_ex_exactly_of_type(pattern, wildcard)) {
 474
 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);
 481                 }
 482                 repl_lst.append(pattern == *this);
 483                 return true;
 484
 485         } else {
 486
 487                 // Expression must be of the same type as the pattern
 488                 if (tinfo() != pattern.bp->tinfo())
 489                         return false;
 490
 491                 // Number of subexpressions must match
 492                 if (nops() != pattern.nops())
 493                         return false;
 494
 495                 // No subexpressions? Then just compare the objects (there can't be
 496                 // wildcards in the pattern)
 497                 if (nops() == 0)
 498                         return is_equal_same_type(*pattern.bp);
 499
 500                 // Check whether attributes that are not subexpressions match
 501                 if (!match_same_type(*pattern.bp))
 502                         return false;
 503
 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))
 507                                 return false;
 508
 509                 // Looks similar enough, match found
 510                 return true;
 511         }
 512 }
 513
 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
 517 {
 518         GINAC_ASSERT(ls.nops() == lr.nops());
 519
 520         if (no_pattern) {
 521                 for (unsigned i=0; i<ls.nops(); i++) {
 522                         if (is_equal(*ls.op(i).bp))
 523                                 return lr.op(i);
 524                 }
 525         } else {
 526                 for (unsigned i=0; i<ls.nops(); i++) {
 527                         lst repl_lst;
 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
 530                 }
 531         }
 532
 533         return *this;
 534 }
 535
 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.
 539  *
 540  *  @param s symbol to differentiate in
 541  *  @param nth order of differentiation
 542  *  @see ex::diff */
 543 ex basic::diff(const symbol & s, unsigned nth) const
 544 {
 545         // trivial: zeroth derivative
 546         if (nth==0)
 547                 return ex(*this);
 548
 549         // evaluate unevaluated *this before differentiating
 550         if (!(flags & status_flags::evaluated))
 551                 return ex(*this).diff(s, nth);
 552
 553         ex ndiff = this->derivative(s);
 554         while (!ndiff.is_zero() &&    // stop differentiating zeros
 555                nth>1) {
 556                 ndiff = ndiff.diff(s);
 557                 --nth;
 558         }
 559         return ndiff;
 560 }
 561
 562 /** Return a vector containing the free indices of an expression. */
 563 exvector basic::get_free_indices(void) const
 564 {
 565         return exvector(); // return an empty exvector
 566 }
 567
 568 ex basic::simplify_ncmul(const exvector & v) const
 569 {
 570         return simplified_ncmul(v);
 571 }
 572
 573 // protected
 574
 575 /** Function object to be applied by basic::derivative(). */
 576 struct derivative_map_function : public map_function {
 577         const symbol &s;
 578         derivative_map_function(const symbol &sym) : s(sym) {}
 579         ex operator()(const ex & e) { return diff(e, s); }
 580 };
 581
 582 /** Default implementation of ex::diff(). It maps the operation on the
 583  *  operands (or returns 0 when the object has no operands).
 584  *
 585  *  @see ex::diff */
 586 ex basic::derivative(const symbol & s) const
 587 {
 588         if (nops() == 0)
 589                 return _ex0();
 590         else {
 591                 derivative_map_function map_derivative(s);
 592                 return map(map_derivative);
 593         }
 594 }
 595
 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
 602 {
 603         return compare_pointers(this, &other);
 604 }
 605
 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
 612 {
 613         return compare_same_type(other)==0;
 614 }
 615
 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.
 624  *
 625  *  @see basic::match */
 626 bool basic::match_same_type(const basic & other) const
 627 {
 628         // The default is to only consider subexpressions, but not any other
 629         // attributes
 630         return true;
 631 }
 632
 633 unsigned basic::return_type(void) const
 634 {
 635         return return_types::commutative;
 636 }
 637
 638 unsigned basic::return_type_tinfo(void) const
 639 {
 640         return tinfo();
 641 }
 642
 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
 650 {
 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();
 655         }
 656
 657         // mask out numeric hashes:
 658         v &= 0x7FFFFFFFU;
 659
 660         // store calculated hash value only if object is already evaluated
 661         if (flags & status_flags::evaluated) {
 662                 setflag(status_flags::hash_calculated);
 663                 hashvalue = v;
 664         }
 665
 666         return v;
 667 }
 668
 669 /** Function object to be applied by basic::expand(). */
 670 struct expand_map_function : public map_function {
 671         unsigned options;
 672         expand_map_function(unsigned o) : options(o) {}
 673         ex operator()(const ex & e) { return expand(e, options); }
 674 };
 675
 676 /** Expand expression, i.e. multiply it out and return the result as a new
 677  *  expression. */
 678 ex basic::expand(unsigned options) const
 679 {
 680         if (nops() == 0)
 681                 return (options == 0) ? setflag(status_flags::expanded) : *this;
 682         else {
 683                 expand_map_function map_expand(options);
 684                 return map(map_expand).bp->setflag(options == 0 ? status_flags::expanded : 0);
 685         }
 686 }
 687
 688
 689 //////////
 690 // non-virtual functions in this class
 691 //////////
 692
 693 // public
 694
 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
 701 {
 702         if (e.info(info_flags::relation_equal)) {
 703                 return subs(lst(e), no_pattern);
 704         }
 705         if (!e.info(info_flags::list)) {
 706                 throw(std::invalid_argument("basic::subs(ex): argument must be a list"));
 707         }
 708         lst ls;
 709         lst lr;
 710         for (unsigned i=0; i<e.nops(); i++) {
 711                 ex r = e.op(i);
 712                 if (!r.info(info_flags::relation_equal)) {
 713                         throw(std::invalid_argument("basic::subs(ex): argument must be a list of equations"));
 714                 }
 715                 ls.append(r.op(0));
 716                 lr.append(r.op(1));
 717         }
 718         return subs(ls, lr, no_pattern);
 719 }
 720
 721 /** Compare objects to establish canonical ordering.
 722  *  All compare functions return: -1 for *this less than other, 0 equal,
 723  *  1 greater. */
 724 int basic::compare(const basic & other) const
 725 {
 726         unsigned hash_this = gethash();
 727         unsigned hash_other = other.gethash();
 728
 729         if (hash_this<hash_other) return -1;
 730         if (hash_this>hash_other) return 1;
 731
 732         unsigned typeid_this = tinfo();
 733         unsigned typeid_other = other.tinfo();
 734
 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;
 742                 return -1;
 743         }
 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;
 751                 return 1;
 752         }
 753
 754         GINAC_ASSERT(typeid(*this)==typeid(other));
 755
 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;
 764 //      }
 765 //      return cmpval;
 766
 767         return compare_same_type(other);
 768 }
 769
 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.
 774  *
 775  *  @see is_equal_same_type */
 776 bool basic::is_equal(const basic & other) const
 777 {
 778         if (this->gethash()!=other.gethash())
 779                 return false;
 780         if (this->tinfo()!=other.tinfo())
 781                 return false;
 782
 783         GINAC_ASSERT(typeid(*this)==typeid(other));
 784
 785         return is_equal_same_type(other);
 786 }
 787
 788 // protected
 789
 790 /** Stop further evaluation.
 791  *
 792  *  @see basic::eval */
 793 const basic & basic::hold(void) const
 794 {
 795         return setflag(status_flags::evaluated);
 796 }
 797
 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
 801 {
 802         if (this->refcount>1)
 803                 throw(std::runtime_error("cannot modify multiply referenced object"));
 804         clearflag(status_flags::hash_calculated);
 805 }
 806
 807 //////////
 808 // global variables
 809 //////////
 810
 811 int max_recursion_level = 1024;
 812
 813 } // namespace GiNaC