ginac/basic.cpp

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