ginac/basic.cpp

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