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 "wildcard.h"
  38 #include "print.h"
  39 #include "archive.h"
  40 #include "utils.h"
  41 #include "debugmsg.h"
  42
  43 namespace GiNaC {
  44
  45 GINAC_IMPLEMENT_REGISTERED_CLASS_NO_CTORS(basic, void)
  46
  47 //////////
  48 // default ctor, dtor, copy ctor assignment operator and helpers
  49 //////////
  50
  51 // public
  52
  53 basic::basic(const basic & other) : tinfo_key(TINFO_basic), flags(0), refcount(0)
  54 {
  55         debugmsg("basic copy ctor", LOGLEVEL_CONSTRUCT);
  56         copy(other);
  57 }
  58
  59 const basic & basic::operator=(const basic & other)
  60 {
  61         debugmsg("basic operator=", LOGLEVEL_ASSIGNMENT);
  62         if (this != &other) {
  63                 destroy(true);
  64                 copy(other);
  65         }
  66         return *this;
  67 }
  68
  69 // protected
  70
  71 // none (all conditionally inlined)
  72
  73 //////////
  74 // other ctors
  75 //////////
  76
  77 // none (all conditionally inlined)
  78
  79 //////////
  80 // archiving
  81 //////////
  82
  83 /** Construct object from archive_node. */
  84 basic::basic(const archive_node &n, const lst &sym_lst) : flags(0), refcount(0)
  85 {
  86         debugmsg("basic ctor from archive_node", LOGLEVEL_CONSTRUCT);
  87
  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);
  92         else
  93                 throw (std::runtime_error("archive node contains no class name"));
  94 }
  95
  96 /** Unarchive the object. */
  97 DEFAULT_UNARCHIVE(basic)
  98
  99 /** Archive the object. */
 100 void basic::archive(archive_node &n) const
 101 {
 102         n.add_string("class", class_name());
 103 }
 104
 105 //////////
 106 // new virtual functions which can be overridden by derived classes
 107 //////////
 108
 109 // public
 110
 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
 116 {
 117         debugmsg("basic print", LOGLEVEL_PRINT);
 118
 119         if (is_of_type(c, print_tree)) {
 120
 121                 c.s << std::string(level, ' ') << class_name()
 122                     << std::hex << ", hash=0x" << hashvalue << ", flags=0x" << flags << std::dec
 123                     << ", nops=" << nops()
 124                     << std::endl;
 125                 for (unsigned i=0; i<nops(); ++i)
 126                         op(i).print(c, level + static_cast<const print_tree &>(c).delta_indent);
 127
 128         } else
 129                 c.s << "[" << class_name() << " object]";
 130 }
 131
 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.
 136  *
 137  *  @see basic::print */
 138 void basic::dbgprint(void) const
 139 {
 140         this->print(std::cerr);
 141         std::cerr << std::endl;
 142 }
 143
 144 /** Little wrapper around printtree to be called within a debugger.
 145  *
 146  *  @see basic::dbgprint
 147  *  @see basic::printtree */
 148 void basic::dbgprinttree(void) const
 149 {
 150         this->print(print_tree(std::cerr));
 151 }
 152
 153 /** Return relative operator precedence (for parenthizing output). */
 154 unsigned basic::precedence(void) const
 155 {
 156         return 70;
 157 }
 158
 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
 163 {
 164         debugmsg("basic duplicate",LOGLEVEL_DUPLICATE);
 165         return new basic(*this);
 166 }
 167
 168 /** Information about the object.
 169  *
 170  *  @see class info_flags */
 171 bool basic::info(unsigned inf) const
 172 {
 173         // all possible properties are false for basic objects
 174         return false;
 175 }
 176
 177 /** Number of operands/members. */
 178 unsigned basic::nops() const
 179 {
 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
 182         // override this.
 183         return 0;
 184 }
 185
 186 /** Return operand/member at position i. */
 187 ex basic::op(int i) const
 188 {
 189         return (const_cast<basic *>(this))->let_op(i);
 190 }
 191
 192 /** Return modifyable operand/member at position i. */
 193 ex & basic::let_op(int i)
 194 {
 195         throw(std::out_of_range("op() out of range"));
 196 }
 197
 198 ex basic::operator[](const ex & index) const
 199 {
 200         if (is_ex_exactly_of_type(index,numeric))
 201                 return op(ex_to<numeric>(index).to_int());
 202
 203         throw(std::invalid_argument("non-numeric indices not supported by this type"));
 204 }
 205
 206 ex basic::operator[](int i) const
 207 {
 208         return op(i);
 209 }
 210
 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
 216 {
 217         lst repl_lst;
 218         if (match(pattern, repl_lst))
 219                 return true;
 220         for (unsigned i=0; i<nops(); i++)
 221                 if (op(i).has(pattern))
 222                         return true;
 223
 224         return false;
 225 }
 226
 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
 230 {
 231         unsigned num = nops();
 232         if (num == 0)
 233                 return *this;
 234
 235         basic *copy = duplicate();
 236         copy->setflag(status_flags::dynallocated);
 237         copy->clearflag(status_flags::hash_calculated | status_flags::expanded);
 238         ex e(*copy);
 239         for (unsigned i=0; i<num; i++)
 240                 e.let_op(i) = f(e.op(i));
 241         return e.eval();
 242 }
 243
 244 /** Return degree of highest power in object s. */
 245 int basic::degree(const ex & s) const
 246 {
 247         return 0;
 248 }
 249
 250 /** Return degree of lowest power in object s. */
 251 int basic::ldegree(const ex & s) const
 252 {
 253         return 0;
 254 }
 255
 256 /** Return coefficient of degree n in object s. */
 257 ex basic::coeff(const ex & s, int n) const
 258 {
 259         return n==0 ? *this : _ex0();
 260 }
 261
 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
 266 {
 267         ex x;
 268         if (is_ex_of_type(s, lst)) {
 269
 270                 // List of objects specified
 271                 if (s.nops() == 0)
 272                         return *this;
 273                 if (s.nops() == 1)
 274                         return collect(s.op(0));
 275
 276                 else if (distributed) {
 277
 278                         // Get lower/upper degree of all symbols in list
 279                         int num = s.nops();
 280                         struct sym_info {
 281                                 ex sym;
 282                                 int ldeg, deg;
 283                                 int cnt;  // current degree, 'counter'
 284                                 ex coeff; // coefficient for degree 'cnt'
 285                         };
 286                         sym_info *si = new sym_info[num];
 287                         ex c = *this;
 288                         for (int i=0; i<num; i++) {
 289                                 si[i].sym = s.op(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);
 293                         }
 294
 295                         while (true) {
 296
 297                                 // Calculate coeff*x1^c1*...*xn^cn
 298                                 ex y = _ex1();
 299                                 for (int i=0; i<num; i++) {
 300                                         int cnt = si[i].cnt;
 301                                         y *= power(si[i].sym, cnt);
 302                                 }
 303                                 x += y * si[num - 1].coeff;
 304
 305                                 // Increment counters
 306                                 int n = num - 1;
 307                                 while (true) {
 308                                         ++si[n].cnt;
 309                                         if (si[n].cnt <= si[n].deg) {
 310                                                 // Update coefficients
 311                                                 ex c;
 312                                                 if (n == 0)
 313                                                         c = *this;
 314                                                 else
 315                                                         c = si[n - 1].coeff;
 316                                                 for (int i=n; i<num; i++)
 317                                                         c = si[i].coeff = c.coeff(si[i].sym, si[i].cnt);
 318                                                 break;
 319                                         }
 320                                         if (n == 0)
 321                                                 goto done;
 322                                         si[n].cnt = si[n].ldeg;
 323                                         n--;
 324                                 }
 325                         }
 326
 327 done:           delete[] si;
 328
 329                 } else {
 330
 331                         // Recursive form
 332                         x = *this;
 333                         for (int n=s.nops()-1; n>=0; n--)
 334                                 x = x.collect(s[n]);
 335                 }
 336
 337         } else {
 338
 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);
 342         }
 343
 344         // correct for lost fractional arguments and return
 345         return x + (*this - x).expand();
 346 }
 347
 348 /** Perform automatic non-interruptive term rewriting rules. */
 349 ex basic::eval(int level) const
 350 {
 351         // There is nothing to do for basic objects:
 352         return this->hold();
 353 }
 354
 355 /** Function object to be applied by basic::evalf(). */
 356 struct evalf_map_function : public map_function {
 357         int level;
 358         evalf_map_function(int l) : level(l) {}
 359         ex operator()(const ex & e) { return evalf(e, level); }
 360 };
 361
 362 /** Evaluate object numerically. */
 363 ex basic::evalf(int level) const
 364 {
 365         if (nops() == 0)
 366                 return *this;
 367         else {
 368                 if (level == 1)
 369                         return *this;
 370                 else if (level == -max_recursion_level)
 371                         throw(std::runtime_error("max recursion level reached"));
 372                 else {
 373                         evalf_map_function map_evalf(level - 1);
 374                         return map(map_evalf);
 375                 }
 376         }
 377 }
 378
 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); }
 382 } map_evalm;
 383
 384 /** Evaluate sums, products and integer powers of matrices. */
 385 ex basic::evalm(void) const
 386 {
 387         if (nops() == 0)
 388                 return *this;
 389         else
 390                 return map(map_evalm);
 391 }
 392
 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
 398 {
 399         // There is nothing to do for basic objects
 400         return i.hold();
 401 }
 402
 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().
 406  *
 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
 412 {
 413         return self + other;
 414 }
 415
 416 /** Multiply an indexed expression with a scalar. This function is used
 417  *  internally by simplify_indexed().
 418  *
 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
 424 {
 425         return self * other;
 426 }
 427
 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().
 434  *
 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
 441 {
 442         // Do nothing
 443         return false;
 444 }
 445
 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
 450 {
 451 /*
 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.
 465 */
 466
 467         if (is_ex_exactly_of_type(pattern, wildcard)) {
 468
 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)));
 475                 }
 476                 repl_lst.append(pattern == *this);
 477                 return true;
 478
 479         } else {
 480
 481                 // Expression must be of the same type as the pattern
 482                 if (tinfo() != ex_to<basic>(pattern).tinfo())
 483                         return false;
 484
 485                 // Number of subexpressions must match
 486                 if (nops() != pattern.nops())
 487                         return false;
 488
 489                 // No subexpressions? Then just compare the objects (there can't be
 490                 // wildcards in the pattern)
 491                 if (nops() == 0)
 492                         return is_equal_same_type(ex_to<basic>(pattern));
 493
 494                 // Check whether attributes that are not subexpressions match
 495                 if (!match_same_type(ex_to<basic>(pattern)))
 496                         return false;
 497
 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))
 501                                 return false;
 502
 503                 // Looks similar enough, match found
 504                 return true;
 505         }
 506 }
 507
 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
 511 {
 512         GINAC_ASSERT(ls.nops() == lr.nops());
 513
 514         if (no_pattern) {
 515                 for (unsigned i=0; i<ls.nops(); i++) {
 516                         if (is_equal(ex_to<basic>(ls.op(i))))
 517                                 return lr.op(i);
 518                 }
 519         } else {
 520                 for (unsigned i=0; i<ls.nops(); i++) {
 521                         lst repl_lst;
 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
 524                 }
 525         }
 526
 527         return *this;
 528 }
 529
 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.
 533  *
 534  *  @param s symbol to differentiate in
 535  *  @param nth order of differentiation
 536  *  @see ex::diff */
 537 ex basic::diff(const symbol & s, unsigned nth) const
 538 {
 539         // trivial: zeroth derivative
 540         if (nth==0)
 541                 return ex(*this);
 542
 543         // evaluate unevaluated *this before differentiating
 544         if (!(flags & status_flags::evaluated))
 545                 return ex(*this).diff(s, nth);
 546
 547         ex ndiff = this->derivative(s);
 548         while (!ndiff.is_zero() &&    // stop differentiating zeros
 549                nth>1) {
 550                 ndiff = ndiff.diff(s);
 551                 --nth;
 552         }
 553         return ndiff;
 554 }
 555
 556 /** Return a vector containing the free indices of an expression. */
 557 exvector basic::get_free_indices(void) const
 558 {
 559         return exvector(); // return an empty exvector
 560 }
 561
 562 ex basic::simplify_ncmul(const exvector & v) const
 563 {
 564         return simplified_ncmul(v);
 565 }
 566
 567 // protected
 568
 569 /** Function object to be applied by basic::derivative(). */
 570 struct derivative_map_function : public map_function {
 571         const symbol &s;
 572         derivative_map_function(const symbol &sym) : s(sym) {}
 573         ex operator()(const ex & e) { return diff(e, s); }
 574 };
 575
 576 /** Default implementation of ex::diff(). It maps the operation on the
 577  *  operands (or returns 0 when the object has no operands).
 578  *
 579  *  @see ex::diff */
 580 ex basic::derivative(const symbol & s) const
 581 {
 582         if (nops() == 0)
 583                 return _ex0();
 584         else {
 585                 derivative_map_function map_derivative(s);
 586                 return map(map_derivative);
 587         }
 588 }
 589
 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
 596 {
 597         return compare_pointers(this, &other);
 598 }
 599
 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
 606 {
 607         return compare_same_type(other)==0;
 608 }
 609
 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.
 618  *
 619  *  @see basic::match */
 620 bool basic::match_same_type(const basic & other) const
 621 {
 622         // The default is to only consider subexpressions, but not any other
 623         // attributes
 624         return true;
 625 }
 626
 627 unsigned basic::return_type(void) const
 628 {
 629         return return_types::commutative;
 630 }
 631
 632 unsigned basic::return_type_tinfo(void) const
 633 {
 634         return tinfo();
 635 }
 636
 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
 644 {
 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();
 649         }
 650
 651         // mask out numeric hashes:
 652         v &= 0x7FFFFFFFU;
 653
 654         // store calculated hash value only if object is already evaluated
 655         if (flags & status_flags::evaluated) {
 656                 setflag(status_flags::hash_calculated);
 657                 hashvalue = v;
 658         }
 659
 660         return v;
 661 }
 662
 663 /** Function object to be applied by basic::expand(). */
 664 struct expand_map_function : public map_function {
 665         unsigned options;
 666         expand_map_function(unsigned o) : options(o) {}
 667         ex operator()(const ex & e) { return expand(e, options); }
 668 };
 669
 670 /** Expand expression, i.e. multiply it out and return the result as a new
 671  *  expression. */
 672 ex basic::expand(unsigned options) const
 673 {
 674         if (nops() == 0)
 675                 return (options == 0) ? setflag(status_flags::expanded) : *this;
 676         else {
 677                 expand_map_function map_expand(options);
 678                 return ex_to<basic>(map(map_expand)).setflag(options == 0 ? status_flags::expanded : 0);
 679         }
 680 }
 681
 682
 683 //////////
 684 // non-virtual functions in this class
 685 //////////
 686
 687 // public
 688
 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
 695 {
 696         if (e.info(info_flags::relation_equal)) {
 697                 return subs(lst(e), no_pattern);
 698         }
 699         if (!e.info(info_flags::list)) {
 700                 throw(std::invalid_argument("basic::subs(ex): argument must be a list"));
 701         }
 702         lst ls;
 703         lst lr;
 704         for (unsigned i=0; i<e.nops(); i++) {
 705                 ex r = e.op(i);
 706                 if (!r.info(info_flags::relation_equal)) {
 707                         throw(std::invalid_argument("basic::subs(ex): argument must be a list of equations"));
 708                 }
 709                 ls.append(r.op(0));
 710                 lr.append(r.op(1));
 711         }
 712         return subs(ls, lr, no_pattern);
 713 }
 714
 715 /** Compare objects to establish canonical ordering.
 716  *  All compare functions return: -1 for *this less than other, 0 equal,
 717  *  1 greater. */
 718 int basic::compare(const basic & other) const
 719 {
 720         unsigned hash_this = gethash();
 721         unsigned hash_other = other.gethash();
 722
 723         if (hash_this<hash_other) return -1;
 724         if (hash_this>hash_other) return 1;
 725
 726         unsigned typeid_this = tinfo();
 727         unsigned typeid_other = other.tinfo();
 728
 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;
 736                 return -1;
 737         }
 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;
 745                 return 1;
 746         }
 747
 748         GINAC_ASSERT(typeid(*this)==typeid(other));
 749
 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;
 758 //      }
 759 //      return cmpval;
 760
 761         return compare_same_type(other);
 762 }
 763
 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.
 768  *
 769  *  @see is_equal_same_type */
 770 bool basic::is_equal(const basic & other) const
 771 {
 772         if (this->gethash()!=other.gethash())
 773                 return false;
 774         if (this->tinfo()!=other.tinfo())
 775                 return false;
 776
 777         GINAC_ASSERT(typeid(*this)==typeid(other));
 778
 779         return is_equal_same_type(other);
 780 }
 781
 782 // protected
 783
 784 /** Stop further evaluation.
 785  *
 786  *  @see basic::eval */
 787 const basic & basic::hold(void) const
 788 {
 789         return setflag(status_flags::evaluated);
 790 }
 791
 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
 795 {
 796         if (this->refcount>1)
 797                 throw(std::runtime_error("cannot modify multiply referenced object"));
 798         clearflag(status_flags::hash_calculated);
 799 }
 800
 801 //////////
 802 // global variables
 803 //////////
 804
 805 int max_recursion_level = 1024;
 806
 807 } // namespace GiNaC