ginac/basic.cpp

   1 /** @file basic.cpp
   2  *
   3  *  Implementation of GiNaC's ABC. */
   4
   5 /*
   6  *  GiNaC Copyright (C) 1999-2003 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 "operators.h"
  38 #include "wildcard.h"
  39 #include "print.h"
  40 #include "archive.h"
  41 #include "utils.h"
  42
  43 namespace GiNaC {
  44
  45 GINAC_IMPLEMENT_REGISTERED_CLASS(basic, void)
  46
  47 //////////
  48 // default constructor, destructor, copy constructor and assignment operator
  49 //////////
  50
  51 // public
  52
  53 /** basic copy constructor: implicitly assumes that the other class is of
  54  *  the exact same type (as it's used by duplicate()), so it can copy the
  55  *  tinfo_key and the hash value. */
  56 basic::basic(const basic & other) : tinfo_key(other.tinfo_key), flags(other.flags & ~status_flags::dynallocated), hashvalue(other.hashvalue), refcount(0)
  57 {
  58         GINAC_ASSERT(typeid(*this) == typeid(other));
  59 }
  60
  61 /** basic assignment operator: the other object might be of a derived class. */
  62 const basic & basic::operator=(const basic & other)
  63 {
  64         unsigned fl = other.flags & ~status_flags::dynallocated;
  65         if (tinfo_key != other.tinfo_key) {
  66                 // The other object is of a derived class, so clear the flags as they
  67                 // might no longer apply (especially hash_calculated). Oh, and don't
  68                 // copy the tinfo_key: it is already set correctly for this object.
  69                 flags = 0;
  70         } else {
  71                 // The objects are of the exact same class, so copy the hash value.
  72                 hashvalue = other.hashvalue;
  73         }
  74         flags = fl;
  75         refcount = 0;
  76         return *this;
  77 }
  78
  79 // protected
  80
  81 // none (all inlined)
  82
  83 //////////
  84 // other constructors
  85 //////////
  86
  87 // none (all inlined)
  88
  89 //////////
  90 // archiving
  91 //////////
  92
  93 /** Construct object from archive_node. */
  94 basic::basic(const archive_node &n, lst &sym_lst) : flags(0), refcount(0)
  95 {
  96         // Reconstruct tinfo_key from class name
  97         std::string class_name;
  98         if (n.find_string("class", class_name))
  99                 tinfo_key = find_tinfo_key(class_name);
 100         else
 101                 throw (std::runtime_error("archive node contains no class name"));
 102 }
 103
 104 /** Unarchive the object. */
 105 DEFAULT_UNARCHIVE(basic)
 106
 107 /** Archive the object. */
 108 void basic::archive(archive_node &n) const
 109 {
 110         n.add_string("class", class_name());
 111 }
 112
 113 //////////
 114 // new virtual functions which can be overridden by derived classes
 115 //////////
 116
 117 // public
 118
 119 /** Output to stream.
 120  *  @param c print context object that describes the output formatting
 121  *  @param level value that is used to identify the precedence or indentation
 122  *               level for placing parentheses and formatting */
 123 void basic::print(const print_context & c, unsigned level) const
 124 {
 125         if (is_a<print_tree>(c)) {
 126
 127                 c.s << std::string(level, ' ') << class_name()
 128                     << std::hex << ", hash=0x" << hashvalue << ", flags=0x" << flags << std::dec
 129                     << ", nops=" << nops()
 130                     << std::endl;
 131                 for (size_t i=0; i<nops(); ++i)
 132                         op(i).print(c, level + static_cast<const print_tree &>(c).delta_indent);
 133
 134         } else
 135                 c.s << "[" << class_name() << " object]";
 136 }
 137
 138 /** Little wrapper around print to be called within a debugger.
 139  *  This is needed because you cannot call foo.print(cout) from within the
 140  *  debugger because it might not know what cout is.  This method can be
 141  *  invoked with no argument and it will simply print to stdout.
 142  *
 143  *  @see basic::print
 144  *  @see basic::dbgprinttree */
 145 void basic::dbgprint() const
 146 {
 147         this->print(std::cerr);
 148         std::cerr << std::endl;
 149 }
 150
 151 /** Little wrapper around printtree to be called within a debugger.
 152  *
 153  *  @see basic::dbgprint */
 154 void basic::dbgprinttree() const
 155 {
 156         this->print(print_tree(std::cerr));
 157 }
 158
 159 /** Return relative operator precedence (for parenthizing output). */
 160 unsigned basic::precedence() const
 161 {
 162         return 70;
 163 }
 164
 165 /** Information about the object.
 166  *
 167  *  @see class info_flags */
 168 bool basic::info(unsigned inf) const
 169 {
 170         // all possible properties are false for basic objects
 171         return false;
 172 }
 173
 174 /** Number of operands/members. */
 175 size_t basic::nops() const
 176 {
 177         // iterating from 0 to nops() on atomic objects should be an empty loop,
 178         // and accessing their elements is a range error.  Container objects should
 179         // override this.
 180         return 0;
 181 }
 182
 183 /** Return operand/member at position i. */
 184 ex basic::op(size_t i) const
 185 {
 186         throw(std::range_error(std::string("basic::op(): ") + class_name() + std::string(" has no operands")));
 187 }
 188
 189 /** Return modifyable operand/member at position i. */
 190 ex & basic::let_op(size_t i)
 191 {
 192         ensure_if_modifiable();
 193         throw(std::range_error(std::string("basic::let_op(): ") + class_name() + std::string(" has no operands")));
 194 }
 195
 196 ex basic::operator[](const ex & index) const
 197 {
 198         if (is_exactly_a<numeric>(index))
 199                 return op(static_cast<size_t>(ex_to<numeric>(index).to_int()));
 200
 201         throw(std::invalid_argument(std::string("non-numeric indices not supported by ") + class_name()));
 202 }
 203
 204 ex basic::operator[](size_t i) const
 205 {
 206         return op(i);
 207 }
 208
 209 ex & basic::operator[](const ex & index)
 210 {
 211         if (is_exactly_a<numeric>(index))
 212                 return let_op(ex_to<numeric>(index).to_int());
 213
 214         throw(std::invalid_argument(std::string("non-numeric indices not supported by ") + class_name()));
 215 }
 216
 217 ex & basic::operator[](size_t i)
 218 {
 219         return let_op(i);
 220 }
 221
 222 /** Test for occurrence of a pattern.  An object 'has' a pattern if it matches
 223  *  the pattern itself or one of the children 'has' it.  As a consequence
 224  *  (according to the definition of children) given e=x+y+z, e.has(x) is true
 225  *  but e.has(x+y) is false. */
 226 bool basic::has(const ex & pattern) const
 227 {
 228         lst repl_lst;
 229         if (match(pattern, repl_lst))
 230                 return true;
 231         for (size_t i=0; i<nops(); i++)
 232                 if (op(i).has(pattern))
 233                         return true;
 234
 235         return false;
 236 }
 237
 238 /** Construct new expression by applying the specified function to all
 239  *  sub-expressions (one level only, not recursively). */
 240 ex basic::map(map_function & f) const
 241 {
 242         size_t num = nops();
 243         if (num == 0)
 244                 return *this;
 245
 246         basic *copy = duplicate();
 247         copy->setflag(status_flags::dynallocated);
 248         copy->clearflag(status_flags::hash_calculated | status_flags::expanded);
 249         for (size_t i=0; i<num; i++)
 250                 copy->let_op(i) = f(copy->op(i));
 251         return *copy;
 252 }
 253
 254 /** Return degree of highest power in object s. */
 255 int basic::degree(const ex & s) const
 256 {
 257         return is_equal(ex_to<basic>(s)) ? 1 : 0;
 258 }
 259
 260 /** Return degree of lowest power in object s. */
 261 int basic::ldegree(const ex & s) const
 262 {
 263         return is_equal(ex_to<basic>(s)) ? 1 : 0;
 264 }
 265
 266 /** Return coefficient of degree n in object s. */
 267 ex basic::coeff(const ex & s, int n) const
 268 {
 269         if (is_equal(ex_to<basic>(s)))
 270                 return n==1 ? _ex1 : _ex0;
 271         else
 272                 return n==0 ? *this : _ex0;
 273 }
 274
 275 /** Sort expanded expression in terms of powers of some object(s).
 276  *  @param s object(s) to sort in
 277  *  @param distributed recursive or distributed form (only used when s is a list) */
 278 ex basic::collect(const ex & s, bool distributed) const
 279 {
 280         ex x;
 281         if (is_a<lst>(s)) {
 282
 283                 // List of objects specified
 284                 if (s.nops() == 0)
 285                         return *this;
 286                 if (s.nops() == 1)
 287                         return collect(s.op(0));
 288
 289                 else if (distributed) {
 290
 291                         // Get lower/upper degree of all symbols in list
 292                         size_t num = s.nops();
 293                         struct sym_info {
 294                                 ex sym;
 295                                 int ldeg, deg;
 296                                 int cnt;  // current degree, 'counter'
 297                                 ex coeff; // coefficient for degree 'cnt'
 298                         };
 299                         sym_info *si = new sym_info[num];
 300                         ex c = *this;
 301                         for (size_t i=0; i<num; i++) {
 302                                 si[i].sym = s.op(i);
 303                                 si[i].ldeg = si[i].cnt = this->ldegree(si[i].sym);
 304                                 si[i].deg = this->degree(si[i].sym);
 305                                 c = si[i].coeff = c.coeff(si[i].sym, si[i].cnt);
 306                         }
 307
 308                         while (true) {
 309
 310                                 // Calculate coeff*x1^c1*...*xn^cn
 311                                 ex y = _ex1;
 312                                 for (size_t i=0; i<num; i++) {
 313                                         int cnt = si[i].cnt;
 314                                         y *= power(si[i].sym, cnt);
 315                                 }
 316                                 x += y * si[num - 1].coeff;
 317
 318                                 // Increment counters
 319                                 size_t n = num - 1;
 320                                 while (true) {
 321                                         ++si[n].cnt;
 322                                         if (si[n].cnt <= si[n].deg) {
 323                                                 // Update coefficients
 324                                                 ex c;
 325                                                 if (n == 0)
 326                                                         c = *this;
 327                                                 else
 328                                                         c = si[n - 1].coeff;
 329                                                 for (size_t i=n; i<num; i++)
 330                                                         c = si[i].coeff = c.coeff(si[i].sym, si[i].cnt);
 331                                                 break;
 332                                         }
 333                                         if (n == 0)
 334                                                 goto done;
 335                                         si[n].cnt = si[n].ldeg;
 336                                         n--;
 337                                 }
 338                         }
 339
 340 done:           delete[] si;
 341
 342                 } else {
 343
 344                         // Recursive form
 345                         x = *this;
 346                         size_t n = s.nops() - 1;
 347                         while (true) {
 348                                 x = x.collect(s[n]);
 349                                 if (n == 0)
 350                                         break;
 351                                 n--;
 352                         }
 353                 }
 354
 355         } else {
 356
 357                 // Only one object specified
 358                 for (int n=this->ldegree(s); n<=this->degree(s); ++n)
 359                         x += this->coeff(s,n)*power(s,n);
 360         }
 361
 362         // correct for lost fractional arguments and return
 363         return x + (*this - x).expand();
 364 }
 365
 366 /** Perform automatic non-interruptive term rewriting rules. */
 367 ex basic::eval(int level) const
 368 {
 369         // There is nothing to do for basic objects:
 370         return hold();
 371 }
 372
 373 /** Function object to be applied by basic::evalf(). */
 374 struct evalf_map_function : public map_function {
 375         int level;
 376         evalf_map_function(int l) : level(l) {}
 377         ex operator()(const ex & e) { return evalf(e, level); }
 378 };
 379
 380 /** Evaluate object numerically. */
 381 ex basic::evalf(int level) const
 382 {
 383         if (nops() == 0)
 384                 return *this;
 385         else {
 386                 if (level == 1)
 387                         return *this;
 388                 else if (level == -max_recursion_level)
 389                         throw(std::runtime_error("max recursion level reached"));
 390                 else {
 391                         evalf_map_function map_evalf(level - 1);
 392                         return map(map_evalf);
 393                 }
 394         }
 395 }
 396
 397 /** Function object to be applied by basic::evalm(). */
 398 struct evalm_map_function : public map_function {
 399         ex operator()(const ex & e) { return evalm(e); }
 400 } map_evalm;
 401
 402 /** Evaluate sums, products and integer powers of matrices. */
 403 ex basic::evalm() const
 404 {
 405         if (nops() == 0)
 406                 return *this;
 407         else
 408                 return map(map_evalm);
 409 }
 410
 411 /** Perform automatic symbolic evaluations on indexed expression that
 412  *  contains this object as the base expression. */
 413 ex basic::eval_indexed(const basic & i) const
 414  // this function can't take a "const ex & i" because that would result
 415  // in an infinite eval() loop
 416 {
 417         // There is nothing to do for basic objects
 418         return i.hold();
 419 }
 420
 421 /** Add two indexed expressions. They are guaranteed to be of class indexed
 422  *  (or a subclass) and their indices are compatible. This function is used
 423  *  internally by simplify_indexed().
 424  *
 425  *  @param self First indexed expression; it's base object is *this
 426  *  @param other Second indexed expression
 427  *  @return sum of self and other
 428  *  @see ex::simplify_indexed() */
 429 ex basic::add_indexed(const ex & self, const ex & other) const
 430 {
 431         return self + other;
 432 }
 433
 434 /** Multiply an indexed expression with a scalar. This function is used
 435  *  internally by simplify_indexed().
 436  *
 437  *  @param self Indexed expression; it's base object is *this
 438  *  @param other Numeric value
 439  *  @return product of self and other
 440  *  @see ex::simplify_indexed() */
 441 ex basic::scalar_mul_indexed(const ex & self, const numeric & other) const
 442 {
 443         return self * other;
 444 }
 445
 446 /** Try to contract two indexed expressions that appear in the same product.
 447  *  If a contraction exists, the function overwrites one or both of the
 448  *  expressions and returns true. Otherwise it returns false. It is
 449  *  guaranteed that both expressions are of class indexed (or a subclass)
 450  *  and that at least one dummy index has been found. This functions is
 451  *  used internally by simplify_indexed().
 452  *
 453  *  @param self Pointer to first indexed expression; it's base object is *this
 454  *  @param other Pointer to second indexed expression
 455  *  @param v The complete vector of factors
 456  *  @return true if the contraction was successful, false otherwise
 457  *  @see ex::simplify_indexed() */
 458 bool basic::contract_with(exvector::iterator self, exvector::iterator other, exvector & v) const
 459 {
 460         // Do nothing
 461         return false;
 462 }
 463
 464 /** Check whether the expression matches a given pattern. For every wildcard
 465  *  object in the pattern, an expression of the form "wildcard == matching_expression"
 466  *  is added to repl_lst. */
 467 bool basic::match(const ex & pattern, lst & repl_lst) const
 468 {
 469 /*
 470         Sweet sweet shapes, sweet sweet shapes,
 471         That's the key thing, right right.
 472         Feed feed face, feed feed shapes,
 473         But who is the king tonight?
 474         Who is the king tonight?
 475         Pattern is the thing, the key thing-a-ling,
 476         But who is the king of Pattern?
 477         But who is the king, the king thing-a-ling,
 478         Who is the king of Pattern?
 479         Bog is the king, the king thing-a-ling,
 480         Bog is the king of Pattern.
 481         Ba bu-bu-bu-bu bu-bu-bu-bu-bu-bu bu-bu
 482         Bog is the king of Pattern.
 483 */
 484
 485         if (is_exactly_a<wildcard>(pattern)) {
 486
 487                 // Wildcard matches anything, but check whether we already have found
 488                 // a match for that wildcard first (if so, it the earlier match must
 489                 // be the same expression)
 490                 for (lst::const_iterator it = repl_lst.begin(); it != repl_lst.end(); ++it) {
 491                         if (it->op(0).is_equal(pattern))
 492                                 return is_equal(ex_to<basic>(it->op(1)));
 493                 }
 494                 repl_lst.append(pattern == *this);
 495                 return true;
 496
 497         } else {
 498
 499                 // Expression must be of the same type as the pattern
 500                 if (tinfo() != ex_to<basic>(pattern).tinfo())
 501                         return false;
 502
 503                 // Number of subexpressions must match
 504                 if (nops() != pattern.nops())
 505                         return false;
 506
 507                 // No subexpressions? Then just compare the objects (there can't be
 508                 // wildcards in the pattern)
 509                 if (nops() == 0)
 510                         return is_equal_same_type(ex_to<basic>(pattern));
 511
 512                 // Check whether attributes that are not subexpressions match
 513                 if (!match_same_type(ex_to<basic>(pattern)))
 514                         return false;
 515
 516                 // Otherwise the subexpressions must match one-to-one
 517                 for (size_t i=0; i<nops(); i++)
 518                         if (!op(i).match(pattern.op(i), repl_lst))
 519                                 return false;
 520
 521                 // Looks similar enough, match found
 522                 return true;
 523         }
 524 }
 525
 526 /** Helper function for subs(). Does not recurse into subexpressions. */
 527 ex basic::subs_one_level(const lst & ls, const lst & lr, unsigned options) const
 528 {
 529         GINAC_ASSERT(ls.nops() == lr.nops());
 530
 531         lst::const_iterator its, itr;
 532
 533         if (options & subs_options::subs_no_pattern) {
 534                 for (its = ls.begin(), itr = lr.begin(); its != ls.end(); ++its, ++itr) {
 535                         if (is_equal(ex_to<basic>(*its)))
 536                                 return *itr;
 537                 }
 538         } else {
 539                 for (its = ls.begin(), itr = lr.begin(); its != ls.end(); ++its, ++itr) {
 540                         lst repl_lst;
 541                         if (match(ex_to<basic>(*its), repl_lst))
 542                                 return itr->subs(repl_lst, options | subs_options::subs_no_pattern); // avoid infinite recursion when re-substituting the wildcards
 543                 }
 544         }
 545
 546         return *this;
 547 }
 548
 549 /** Substitute a set of objects by arbitrary expressions. The ex returned
 550  *  will already be evaluated. */
 551 ex basic::subs(const lst & ls, const lst & lr, unsigned options) const
 552 {
 553         size_t num = nops();
 554         if (num) {
 555
 556                 // Substitute in subexpressions
 557                 for (size_t i=0; i<num; i++) {
 558                         const ex & orig_op = op(i);
 559                         const ex & subsed_op = orig_op.subs(ls, lr, options);
 560                         if (!are_ex_trivially_equal(orig_op, subsed_op)) {
 561
 562                                 // Something changed, clone the object
 563                                 basic *copy = duplicate();
 564                                 copy->setflag(status_flags::dynallocated);
 565                                 copy->clearflag(status_flags::hash_calculated | status_flags::expanded);
 566
 567                                 // Substitute the changed operand
 568                                 copy->let_op(i++) = subsed_op;
 569
 570                                 // Substitute the other operands
 571                                 for (; i<num; i++)
 572                                         copy->let_op(i) = op(i).subs(ls, lr, options);
 573
 574                                 // Perform substitutions on the new object as a whole
 575                                 return copy->subs_one_level(ls, lr, options);
 576                         }
 577                 }
 578         }
 579
 580         // Nothing changed or no subexpressions
 581         return subs_one_level(ls, lr, options);
 582 }
 583
 584 /** Default interface of nth derivative ex::diff(s, n).  It should be called
 585  *  instead of ::derivative(s) for first derivatives and for nth derivatives it
 586  *  just recurses down.
 587  *
 588  *  @param s symbol to differentiate in
 589  *  @param nth order of differentiation
 590  *  @see ex::diff */
 591 ex basic::diff(const symbol & s, unsigned nth) const
 592 {
 593         // trivial: zeroth derivative
 594         if (nth==0)
 595                 return ex(*this);
 596
 597         // evaluate unevaluated *this before differentiating
 598         if (!(flags & status_flags::evaluated))
 599                 return ex(*this).diff(s, nth);
 600
 601         ex ndiff = this->derivative(s);
 602         while (!ndiff.is_zero() &&    // stop differentiating zeros
 603                nth>1) {
 604                 ndiff = ndiff.diff(s);
 605                 --nth;
 606         }
 607         return ndiff;
 608 }
 609
 610 /** Return a vector containing the free indices of an expression. */
 611 exvector basic::get_free_indices() const
 612 {
 613         return exvector(); // return an empty exvector
 614 }
 615
 616 ex basic::eval_ncmul(const exvector & v) const
 617 {
 618         return hold_ncmul(v);
 619 }
 620
 621 // protected
 622
 623 /** Function object to be applied by basic::derivative(). */
 624 struct derivative_map_function : public map_function {
 625         const symbol &s;
 626         derivative_map_function(const symbol &sym) : s(sym) {}
 627         ex operator()(const ex & e) { return diff(e, s); }
 628 };
 629
 630 /** Default implementation of ex::diff(). It maps the operation on the
 631  *  operands (or returns 0 when the object has no operands).
 632  *
 633  *  @see ex::diff */
 634 ex basic::derivative(const symbol & s) const
 635 {
 636         if (nops() == 0)
 637                 return _ex0;
 638         else {
 639                 derivative_map_function map_derivative(s);
 640                 return map(map_derivative);
 641         }
 642 }
 643
 644 /** Returns order relation between two objects of same type.  This needs to be
 645  *  implemented by each class. It may never return anything else than 0,
 646  *  signalling equality, or +1 and -1 signalling inequality and determining
 647  *  the canonical ordering.  (Perl hackers will wonder why C++ doesn't feature
 648  *  the spaceship operator <=> for denoting just this.) */
 649 int basic::compare_same_type(const basic & other) const
 650 {
 651         return compare_pointers(this, &other);
 652 }
 653
 654 /** Returns true if two objects of same type are equal.  Normally needs
 655  *  not be reimplemented as long as it wasn't overwritten by some parent
 656  *  class, since it just calls compare_same_type().  The reason why this
 657  *  function exists is that sometimes it is easier to determine equality
 658  *  than an order relation and then it can be overridden. */
 659 bool basic::is_equal_same_type(const basic & other) const
 660 {
 661         return compare_same_type(other)==0;
 662 }
 663
 664 /** Returns true if the attributes of two objects are similar enough for
 665  *  a match. This function must not match subexpressions (this is already
 666  *  done by basic::match()). Only attributes not accessible by op() should
 667  *  be compared. This is also the reason why this function doesn't take the
 668  *  wildcard replacement list from match() as an argument: only subexpressions
 669  *  are subject to wildcard matches. Also, this function only needs to be
 670  *  implemented for container classes because is_equal_same_type() is
 671  *  automatically used instead of match_same_type() if nops() == 0.
 672  *
 673  *  @see basic::match */
 674 bool basic::match_same_type(const basic & other) const
 675 {
 676         // The default is to only consider subexpressions, but not any other
 677         // attributes
 678         return true;
 679 }
 680
 681 unsigned basic::return_type() const
 682 {
 683         return return_types::commutative;
 684 }
 685
 686 unsigned basic::return_type_tinfo() const
 687 {
 688         return tinfo();
 689 }
 690
 691 /** Compute the hash value of an object and if it makes sense to store it in
 692  *  the objects status_flags, do so.  The method inherited from class basic
 693  *  computes a hash value based on the type and hash values of possible
 694  *  members.  For this reason it is well suited for container classes but
 695  *  atomic classes should override this implementation because otherwise they
 696  *  would all end up with the same hashvalue. */
 697 unsigned basic::calchash() const
 698 {
 699         unsigned v = golden_ratio_hash(tinfo());
 700         for (size_t i=0; i<nops(); i++) {
 701                 v = rotate_left(v);
 702                 v ^= (const_cast<basic *>(this))->op(i).gethash();
 703         }
 704
 705         // store calculated hash value only if object is already evaluated
 706         if (flags & status_flags::evaluated) {
 707                 setflag(status_flags::hash_calculated);
 708                 hashvalue = v;
 709         }
 710
 711         return v;
 712 }
 713
 714 /** Function object to be applied by basic::expand(). */
 715 struct expand_map_function : public map_function {
 716         unsigned options;
 717         expand_map_function(unsigned o) : options(o) {}
 718         ex operator()(const ex & e) { return expand(e, options); }
 719 };
 720
 721 /** Expand expression, i.e. multiply it out and return the result as a new
 722  *  expression. */
 723 ex basic::expand(unsigned options) const
 724 {
 725         if (nops() == 0)
 726                 return (options == 0) ? setflag(status_flags::expanded) : *this;
 727         else {
 728                 expand_map_function map_expand(options);
 729                 return ex_to<basic>(map(map_expand)).setflag(options == 0 ? status_flags::expanded : 0);
 730         }
 731 }
 732
 733
 734 //////////
 735 // non-virtual functions in this class
 736 //////////
 737
 738 // public
 739
 740 /** Substitute objects in an expression (syntactic substitution) and return
 741  *  the result as a new expression.  There are two valid types of
 742  *  replacement arguments: 1) a relational like object==ex and 2) a list of
 743  *  relationals lst(object1==ex1,object2==ex2,...), which is converted to
 744  *  subs(lst(object1,object2,...),lst(ex1,ex2,...)). */
 745 ex basic::subs(const ex & e, unsigned options) const
 746 {
 747         if (e.info(info_flags::relation_equal)) {
 748                 return subs(lst(e), options);
 749         }
 750         if (!e.info(info_flags::list)) {
 751                 throw(std::invalid_argument("basic::subs(ex): argument must be a list"));
 752         }
 753
 754         // Split list into two
 755         lst ls;
 756         lst lr;
 757         GINAC_ASSERT(is_a<lst>(e));
 758         for (lst::const_iterator it = ex_to<lst>(e).begin(); it != ex_to<lst>(e).end(); ++it) {
 759                 ex r = *it;
 760                 if (!r.info(info_flags::relation_equal)) {
 761                         throw(std::invalid_argument("basic::subs(ex): argument must be a list of equations"));
 762                 }
 763                 ls.append(r.op(0));
 764                 lr.append(r.op(1));
 765         }
 766         return subs(ls, lr, options);
 767 }
 768
 769 /** Compare objects syntactically to establish canonical ordering.
 770  *  All compare functions return: -1 for *this less than other, 0 equal,
 771  *  1 greater. */
 772 int basic::compare(const basic & other) const
 773 {
 774         const unsigned hash_this = gethash();
 775         const unsigned hash_other = other.gethash();
 776         if (hash_this<hash_other) return -1;
 777         if (hash_this>hash_other) return 1;
 778
 779         const unsigned typeid_this = tinfo();
 780         const unsigned typeid_other = other.tinfo();
 781         if (typeid_this==typeid_other) {
 782                 GINAC_ASSERT(typeid(*this)==typeid(other));
 783 //              int cmpval = compare_same_type(other);
 784 //              if (cmpval!=0) {
 785 //                      std::cout << "hash collision, same type: "
 786 //                                << *this << " and " << other << std::endl;
 787 //                      this->print(print_tree(std::cout));
 788 //                      std::cout << " and ";
 789 //                      other.print(print_tree(std::cout));
 790 //                      std::cout << std::endl;
 791 //              }
 792 //              return cmpval;
 793                 return compare_same_type(other);
 794         } else {
 795 //              std::cout << "hash collision, different types: "
 796 //                        << *this << " and " << other << std::endl;
 797 //              this->print(print_tree(std::cout));
 798 //              std::cout << " and ";
 799 //              other.print(print_tree(std::cout));
 800 //              std::cout << std::endl;
 801                 return (typeid_this<typeid_other ? -1 : 1);
 802         }
 803 }
 804
 805 /** Test for syntactic equality.
 806  *  This is only a quick test, meaning objects should be in the same domain.
 807  *  You might have to .expand(), .normal() objects first, depending on the
 808  *  domain of your computation, to get a more reliable answer.
 809  *
 810  *  @see is_equal_same_type */
 811 bool basic::is_equal(const basic & other) const
 812 {
 813         if (this->gethash()!=other.gethash())
 814                 return false;
 815         if (this->tinfo()!=other.tinfo())
 816                 return false;
 817
 818         GINAC_ASSERT(typeid(*this)==typeid(other));
 819
 820         return is_equal_same_type(other);
 821 }
 822
 823 // protected
 824
 825 /** Stop further evaluation.
 826  *
 827  *  @see basic::eval */
 828 const basic & basic::hold() const
 829 {
 830         return setflag(status_flags::evaluated);
 831 }
 832
 833 /** Ensure the object may be modified without hurting others, throws if this
 834  *  is not the case. */
 835 void basic::ensure_if_modifiable() const
 836 {
 837         if (refcount > 1)
 838                 throw(std::runtime_error("cannot modify multiply referenced object"));
 839         clearflag(status_flags::hash_calculated | status_flags::evaluated);
 840 }
 841
 842 //////////
 843 // global variables
 844 //////////
 845
 846 int max_recursion_level = 1024;
 847
 848 } // namespace GiNaC