ginac/expairseq.cpp

   1 /** @file expairseq.cpp
   2  *
   3  *  Implementation of sequences of expression pairs. */
   4
   5 /*
   6  *  GiNaC Copyright (C) 1999-2001 Johannes Gutenberg University Mainz, Germany
   7  *
   8  *  This program is free software; you can redistribute it and/or modify
   9  *  it under the terms of the GNU General Public License as published by
  10  *  the Free Software Foundation; either version 2 of the License, or
  11  *  (at your option) any later version.
  12  *
  13  *  This program is distributed in the hope that it will be useful,
  14  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
  15  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  16  *  GNU General Public License for more details.
  17  *
  18  *  You should have received a copy of the GNU General Public License
  19  *  along with this program; if not, write to the Free Software
  20  *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
  21  */
  22
  23 #include <algorithm>
  24 #include <string>
  25 #include <stdexcept>
  26
  27 #include "expairseq.h"
  28 #include "lst.h"
  29 #include "mul.h"
  30 #include "power.h"
  31 #include "relational.h"
  32 #include "wildcard.h"
  33 #include "print.h"
  34 #include "archive.h"
  35 #include "utils.h"
  36
  37 #if EXPAIRSEQ_USE_HASHTAB
  38 #include <cmath>
  39 #endif // EXPAIRSEQ_USE_HASHTAB
  40
  41 namespace GiNaC {
  42
  43
  44 GINAC_IMPLEMENT_REGISTERED_CLASS_NO_CTORS(expairseq, basic)
  45
  46 //////////
  47 // helper classes
  48 //////////
  49
  50 class epp_is_less
  51 {
  52 public:
  53         bool operator()(const epp &lh, const epp &rh) const
  54         {
  55                 return (*lh).is_less(*rh);
  56         }
  57 };
  58
  59 //////////
  60 // default ctor, dtor, copy ctor, assignment operator and helpers
  61 //////////
  62
  63 // public
  64
  65 expairseq::expairseq(const expairseq &other)
  66 {
  67         copy(other);
  68 }
  69
  70 const expairseq &expairseq::operator=(const expairseq &other)
  71 {
  72         if (this != &other) {
  73                 destroy(true);
  74                 copy(other);
  75         }
  76         return *this;
  77 }
  78
  79 // protected
  80
  81 /** For use by copy ctor and assignment operator. */
  82 void expairseq::copy(const expairseq &other)
  83 {
  84         inherited::copy(other);
  85         seq = other.seq;
  86         overall_coeff = other.overall_coeff;
  87 #if EXPAIRSEQ_USE_HASHTAB
  88         // copy hashtab
  89         hashtabsize = other.hashtabsize;
  90         if (hashtabsize!=0) {
  91                 hashmask = other.hashmask;
  92                 hashtab.resize(hashtabsize);
  93                 epvector::const_iterator osb = other.seq.begin();
  94                 for (unsigned i=0; i<hashtabsize; ++i) {
  95                         hashtab[i].clear();
  96                         for (epplist::const_iterator cit=other.hashtab[i].begin();
  97                              cit!=other.hashtab[i].end(); ++cit) {
  98                                 hashtab[i].push_back(seq.begin()+((*cit)-osb));
  99                         }
 100                 }
 101         } else {
 102                 hashtab.clear();
 103         }
 104 #endif // EXPAIRSEQ_USE_HASHTAB
 105 }
 106
 107 DEFAULT_DESTROY(expairseq)
 108
 109 //////////
 110 // other ctors
 111 //////////
 112
 113 expairseq::expairseq(const ex &lh, const ex &rh) : inherited(TINFO_expairseq)
 114 {
 115         construct_from_2_ex(lh,rh);
 116         GINAC_ASSERT(is_canonical());
 117 }
 118
 119 expairseq::expairseq(const exvector &v) : inherited(TINFO_expairseq)
 120 {
 121         construct_from_exvector(v);
 122         GINAC_ASSERT(is_canonical());
 123 }
 124
 125 expairseq::expairseq(const epvector &v, const ex &oc)
 126   : inherited(TINFO_expairseq), overall_coeff(oc)
 127 {
 128         construct_from_epvector(v);
 129         GINAC_ASSERT(is_canonical());
 130 }
 131
 132 expairseq::expairseq(epvector *vp, const ex &oc)
 133   : inherited(TINFO_expairseq), overall_coeff(oc)
 134 {
 135         GINAC_ASSERT(vp!=0);
 136         construct_from_epvector(*vp);
 137         delete vp;
 138         GINAC_ASSERT(is_canonical());
 139 }
 140
 141 //////////
 142 // archiving
 143 //////////
 144
 145 expairseq::expairseq(const archive_node &n, const lst &sym_lst) : inherited(n, sym_lst)
 146 #if EXPAIRSEQ_USE_HASHTAB
 147         , hashtabsize(0)
 148 #endif
 149 {
 150         for (unsigned int i=0; true; i++) {
 151                 ex rest;
 152                 ex coeff;
 153                 if (n.find_ex("rest", rest, sym_lst, i) && n.find_ex("coeff", coeff, sym_lst, i))
 154                         seq.push_back(expair(rest, coeff));
 155                 else
 156                         break;
 157         }
 158         n.find_ex("overall_coeff", overall_coeff, sym_lst);
 159 }
 160
 161 void expairseq::archive(archive_node &n) const
 162 {
 163         inherited::archive(n);
 164         epvector::const_iterator i = seq.begin(), iend = seq.end();
 165         while (i != iend) {
 166                 n.add_ex("rest", i->rest);
 167                 n.add_ex("coeff", i->coeff);
 168                 ++i;
 169         }
 170         n.add_ex("overall_coeff", overall_coeff);
 171 }
 172
 173 DEFAULT_UNARCHIVE(expairseq)
 174
 175 //////////
 176 // functions overriding virtual functions from base classes
 177 //////////
 178
 179 // public
 180
 181 basic *expairseq::duplicate() const
 182 {
 183         return new expairseq(*this);
 184 }
 185
 186 void expairseq::print(const print_context &c, unsigned level) const
 187 {
 188         if (is_of_type(c, print_tree)) {
 189
 190                 unsigned delta_indent = static_cast<const print_tree &>(c).delta_indent;
 191
 192                 c.s << std::string(level, ' ') << class_name()
 193                     << std::hex << ", hash=0x" << hashvalue << ", flags=0x" << flags << std::dec
 194                     << ", nops=" << nops()
 195                     << std::endl;
 196                 unsigned num = seq.size();
 197                 for (unsigned i=0; i<num; ++i) {
 198                         seq[i].rest.print(c, level + delta_indent);
 199                         seq[i].coeff.print(c, level + delta_indent);
 200                         if (i != num - 1)
 201                                 c.s << std::string(level + delta_indent, ' ') << "-----" << std::endl;
 202                 }
 203                 if (!overall_coeff.is_equal(default_overall_coeff())) {
 204                         c.s << std::string(level + delta_indent, ' ') << "-----" << std::endl
 205                             << std::string(level + delta_indent, ' ') << "overall_coeff" << std::endl;
 206                         overall_coeff.print(c, level + delta_indent);
 207                 }
 208                 c.s << std::string(level + delta_indent,' ') << "=====" << std::endl;
 209 #if EXPAIRSEQ_USE_HASHTAB
 210                 c.s << std::string(level + delta_indent,' ')
 211                     << "hashtab size " << hashtabsize << std::endl;
 212                 if (hashtabsize == 0) return;
 213 #define MAXCOUNT 5
 214                 unsigned count[MAXCOUNT+1];
 215                 for (int i=0; i<MAXCOUNT+1; ++i)
 216                         count[i] = 0;
 217                 unsigned this_bin_fill;
 218                 unsigned cum_fill_sq = 0;
 219                 unsigned cum_fill = 0;
 220                 for (unsigned i=0; i<hashtabsize; ++i) {
 221                         this_bin_fill = 0;
 222                         if (hashtab[i].size() > 0) {
 223                                 c.s << std::string(level + delta_indent, ' ')
 224                                     << "bin " << i << " with entries ";
 225                                 for (epplist::const_iterator it=hashtab[i].begin();
 226                                      it!=hashtab[i].end(); ++it) {
 227                                         c.s << *it-seq.begin() << " ";
 228                                         ++this_bin_fill;
 229                                 }
 230                                 c.s << std::endl;
 231                                 cum_fill += this_bin_fill;
 232                                 cum_fill_sq += this_bin_fill*this_bin_fill;
 233                         }
 234                         if (this_bin_fill<MAXCOUNT)
 235                                 ++count[this_bin_fill];
 236                         else
 237                                 ++count[MAXCOUNT];
 238                 }
 239                 unsigned fact = 1;
 240                 double cum_prob = 0;
 241                 double lambda = (1.0*seq.size()) / hashtabsize;
 242                 for (int k=0; k<MAXCOUNT; ++k) {
 243                         if (k>0)
 244                                 fact *= k;
 245                         double prob = std::pow(lambda,k)/fact * std::exp(-lambda);
 246                         cum_prob += prob;
 247                         c.s << std::string(level + delta_indent, ' ') << "bins with " << k << " entries: "
 248                             << int(1000.0*count[k]/hashtabsize)/10.0 << "% (expected: "
 249                             << int(prob*1000)/10.0 << ")" << std::endl;
 250                 }
 251                 c.s << std::string(level + delta_indent, ' ') << "bins with more entries: "
 252                     << int(1000.0*count[MAXCOUNT]/hashtabsize)/10.0 << "% (expected: "
 253                     << int((1-cum_prob)*1000)/10.0 << ")" << std::endl;
 254
 255                 c.s << std::string(level + delta_indent, ' ') << "variance: "
 256                     << 1.0/hashtabsize*cum_fill_sq-(1.0/hashtabsize*cum_fill)*(1.0/hashtabsize*cum_fill)
 257                     << std::endl;
 258                 c.s << std::string(level + delta_indent, ' ') << "average fill: "
 259                     << (1.0*cum_fill)/hashtabsize
 260                     << " (should be equal to " << (1.0*seq.size())/hashtabsize << ")" << std::endl;
 261 #endif // EXPAIRSEQ_USE_HASHTAB
 262
 263         } else {
 264                 c.s << "[[";
 265                 printseq(c, ',', precedence(), level);
 266                 c.s << "]]";
 267         }
 268 }
 269
 270 bool expairseq::info(unsigned inf) const
 271 {
 272         return inherited::info(inf);
 273 }
 274
 275 unsigned expairseq::nops() const
 276 {
 277         if (overall_coeff.is_equal(default_overall_coeff()))
 278                 return seq.size();
 279         else
 280                 return seq.size()+1;
 281 }
 282
 283 ex expairseq::op(int i) const
 284 {
 285         if (unsigned(i)<seq.size())
 286                 return recombine_pair_to_ex(seq[i]);
 287         GINAC_ASSERT(!overall_coeff.is_equal(default_overall_coeff()));
 288         return overall_coeff;
 289 }
 290
 291 ex &expairseq::let_op(int i)
 292 {
 293         throw(std::logic_error("let_op not defined for expairseq and derived classes (add,mul,...)"));
 294 }
 295
 296 ex expairseq::map(map_function &f) const
 297 {
 298         epvector *v = new epvector;
 299         v->reserve(seq.size());
 300
 301         epvector::const_iterator cit = seq.begin(), last = seq.end();
 302         while (cit != last) {
 303                 v->push_back(split_ex_to_pair(f(recombine_pair_to_ex(*cit))));
 304                 ++cit;
 305         }
 306
 307         if (overall_coeff.is_equal(default_overall_coeff()))
 308                 return thisexpairseq(v, default_overall_coeff());
 309         else
 310                 return thisexpairseq(v, f(overall_coeff));
 311 }
 312
 313 /** Perform coefficient-wise automatic term rewriting rules in this class. */
 314 ex expairseq::eval(int level) const
 315 {
 316         if ((level==1) && (flags &status_flags::evaluated))
 317                 return *this;
 318
 319         epvector *vp = evalchildren(level);
 320         if (vp==0)
 321                 return this->hold();
 322
 323         return (new expairseq(vp,overall_coeff))->setflag(status_flags::dynallocated | status_flags::evaluated);
 324 }
 325
 326 bool expairseq::match(const ex & pattern, lst & repl_lst) const
 327 {
 328         // This differs from basic::match() because we want "a+b+c+d" to
 329         // match "d+*+b" with "*" being "a+c", and we want to honor commutativity
 330
 331         if (this->tinfo() == ex_to<basic>(pattern).tinfo()) {
 332
 333                 // Check whether global wildcard (one that matches the "rest of the
 334                 // expression", like "*" above) is present
 335                 bool has_global_wildcard = false;
 336                 ex global_wildcard;
 337                 for (unsigned int i=0; i<pattern.nops(); i++) {
 338                         if (is_ex_exactly_of_type(pattern.op(i), wildcard)) {
 339                                 has_global_wildcard = true;
 340                                 global_wildcard = pattern.op(i);
 341                                 break;
 342                         }
 343                 }
 344
 345                 // Unfortunately, this is an O(N^2) operation because we can't
 346                 // sort the pattern in a useful way...
 347
 348                 // Chop into terms
 349                 exvector ops;
 350                 ops.reserve(nops());
 351                 for (unsigned i=0; i<nops(); i++)
 352                         ops.push_back(op(i));
 353
 354                 // Now, for every term of the pattern, look for a matching term in
 355                 // the expression and remove the match
 356                 for (unsigned i=0; i<pattern.nops(); i++) {
 357                         ex p = pattern.op(i);
 358                         if (has_global_wildcard && p.is_equal(global_wildcard))
 359                                 continue;
 360                         exvector::iterator it = ops.begin(), itend = ops.end();
 361                         while (it != itend) {
 362                                 if (it->match(p, repl_lst)) {
 363                                         ops.erase(it);
 364                                         goto found;
 365                                 }
 366                                 ++it;
 367                         }
 368                         return false; // no match found
 369 found:          ;
 370                 }
 371
 372                 if (has_global_wildcard) {
 373
 374                         // Assign all the remaining terms to the global wildcard (unless
 375                         // it has already been matched before, in which case the matches
 376                         // must be equal)
 377                         unsigned num = ops.size();
 378                         epvector *vp = new epvector();
 379                         vp->reserve(num);
 380                         for (unsigned i=0; i<num; i++)
 381                                 vp->push_back(split_ex_to_pair(ops[i]));
 382                         ex rest = thisexpairseq(vp, default_overall_coeff());
 383                         for (unsigned i=0; i<repl_lst.nops(); i++) {
 384                                 if (repl_lst.op(i).op(0).is_equal(global_wildcard))
 385                                         return rest.is_equal(repl_lst.op(i).op(1));
 386                         }
 387                         repl_lst.append(global_wildcard == rest);
 388                         return true;
 389
 390                 } else {
 391
 392                         // No global wildcard, then the match fails if there are any
 393                         // unmatched terms left
 394                         return ops.empty();
 395                 }
 396         }
 397         return inherited::match(pattern, repl_lst);
 398 }
 399
 400 ex expairseq::subs(const lst &ls, const lst &lr, bool no_pattern) const
 401 {
 402         epvector *vp = subschildren(ls, lr, no_pattern);
 403         if (vp)
 404                 return ex_to<basic>(thisexpairseq(vp, overall_coeff)).basic::subs(ls, lr, no_pattern);
 405         else
 406                 return basic::subs(ls, lr, no_pattern);
 407 }
 408
 409 // protected
 410
 411 int expairseq::compare_same_type(const basic &other) const
 412 {
 413         GINAC_ASSERT(is_a<expairseq>(other));
 414         const expairseq &o = static_cast<const expairseq &>(other);
 415
 416         int cmpval;
 417
 418         // compare number of elements
 419         if (seq.size() != o.seq.size())
 420                 return (seq.size()<o.seq.size()) ? -1 : 1;
 421
 422         // compare overall_coeff
 423         cmpval = overall_coeff.compare(o.overall_coeff);
 424         if (cmpval!=0)
 425                 return cmpval;
 426
 427 #if EXPAIRSEQ_USE_HASHTAB
 428         GINAC_ASSERT(hashtabsize==o.hashtabsize);
 429         if (hashtabsize==0) {
 430 #endif // EXPAIRSEQ_USE_HASHTAB
 431                 epvector::const_iterator cit1 = seq.begin();
 432                 epvector::const_iterator cit2 = o.seq.begin();
 433                 epvector::const_iterator last1 = seq.end();
 434                 epvector::const_iterator last2 = o.seq.end();
 435
 436                 for (; (cit1!=last1)&&(cit2!=last2); ++cit1, ++cit2) {
 437                         cmpval = (*cit1).compare(*cit2);
 438                         if (cmpval!=0) return cmpval;
 439                 }
 440
 441                 GINAC_ASSERT(cit1==last1);
 442                 GINAC_ASSERT(cit2==last2);
 443
 444                 return 0;
 445 #if EXPAIRSEQ_USE_HASHTAB
 446         }
 447
 448         // compare number of elements in each hashtab entry
 449         for (unsigned i=0; i<hashtabsize; ++i) {
 450                 unsigned cursize=hashtab[i].size();
 451                 if (cursize != o.hashtab[i].size())
 452                         return (cursize < o.hashtab[i].size()) ? -1 : 1;
 453         }
 454
 455         // compare individual (sorted) hashtab entries
 456         for (unsigned i=0; i<hashtabsize; ++i) {
 457                 unsigned sz = hashtab[i].size();
 458                 if (sz>0) {
 459                         const epplist &eppl1 = hashtab[i];
 460                         const epplist &eppl2 = o.hashtab[i];
 461                         epplist::const_iterator it1 = eppl1.begin();
 462                         epplist::const_iterator it2 = eppl2.begin();
 463                         while (it1!=eppl1.end()) {
 464                                 cmpval = (*(*it1)).compare(*(*it2));
 465                                 if (cmpval!=0)
 466                                         return cmpval;
 467                                 ++it1;
 468                                 ++it2;
 469                         }
 470                 }
 471         }
 472
 473         return 0; // equal
 474 #endif // EXPAIRSEQ_USE_HASHTAB
 475 }
 476
 477 bool expairseq::is_equal_same_type(const basic &other) const
 478 {
 479         const expairseq &o = static_cast<const expairseq &>(other);
 480
 481         // compare number of elements
 482         if (seq.size()!=o.seq.size())
 483                 return false;
 484
 485         // compare overall_coeff
 486         if (!overall_coeff.is_equal(o.overall_coeff))
 487                 return false;
 488
 489 #if EXPAIRSEQ_USE_HASHTAB
 490         // compare number of elements in each hashtab entry
 491         if (hashtabsize!=o.hashtabsize) {
 492                 std::cout << "this:" << std::endl;
 493                 print(print_tree(std::cout));
 494                 std::cout << "other:" << std::endl;
 495                 other.print(print_tree(std::cout));
 496         }
 497
 498         GINAC_ASSERT(hashtabsize==o.hashtabsize);
 499
 500         if (hashtabsize==0) {
 501 #endif // EXPAIRSEQ_USE_HASHTAB
 502                 epvector::const_iterator cit1 = seq.begin();
 503                 epvector::const_iterator cit2 = o.seq.begin();
 504                 epvector::const_iterator last1 = seq.end();
 505
 506                 while (cit1!=last1) {
 507                         if (!(*cit1).is_equal(*cit2)) return false;
 508                         ++cit1;
 509                         ++cit2;
 510                 }
 511
 512                 return true;
 513 #if EXPAIRSEQ_USE_HASHTAB
 514         }
 515
 516         for (unsigned i=0; i<hashtabsize; ++i) {
 517                 if (hashtab[i].size() != o.hashtab[i].size())
 518                         return false;
 519         }
 520
 521         // compare individual sorted hashtab entries
 522         for (unsigned i=0; i<hashtabsize; ++i) {
 523                 unsigned sz = hashtab[i].size();
 524                 if (sz>0) {
 525                         const epplist &eppl1 = hashtab[i];
 526                         const epplist &eppl2 = o.hashtab[i];
 527                         epplist::const_iterator it1 = eppl1.begin();
 528                         epplist::const_iterator it2 = eppl2.begin();
 529                         while (it1!=eppl1.end()) {
 530                                 if (!(*(*it1)).is_equal(*(*it2))) return false;
 531                                 ++it1;
 532                                 ++it2;
 533                         }
 534                 }
 535         }
 536
 537         return true;
 538 #endif // EXPAIRSEQ_USE_HASHTAB
 539 }
 540
 541 unsigned expairseq::return_type(void) const
 542 {
 543         return return_types::noncommutative_composite;
 544 }
 545
 546 unsigned expairseq::calchash(void) const
 547 {
 548         unsigned v = golden_ratio_hash(this->tinfo());
 549         epvector::const_iterator i = seq.begin(), end = seq.end();
 550         while (i != end) {
 551 #if !EXPAIRSEQ_USE_HASHTAB
 552                 v = rotate_left_31(v); // rotation would spoil commutativity
 553 #endif // EXPAIRSEQ_USE_HASHTAB
 554                 v ^= i->rest.gethash();
 555 #if !EXPAIRSEQ_USE_HASHTAB
 556                 v = rotate_left_31(v);
 557                 v ^= i->coeff.gethash();
 558 #endif // EXPAIRSEQ_USE_HASHTAB
 559                 ++i;
 560         }
 561
 562         v ^= overall_coeff.gethash();
 563         v &= 0x7FFFFFFFU;
 564
 565         // store calculated hash value only if object is already evaluated
 566         if (flags &status_flags::evaluated) {
 567                 setflag(status_flags::hash_calculated);
 568                 hashvalue = v;
 569         }
 570
 571         return v;
 572 }
 573
 574 ex expairseq::expand(unsigned options) const
 575 {
 576         epvector *vp = expandchildren(options);
 577         if (vp == NULL) {
 578                 // The terms have not changed, so it is safe to declare this expanded
 579                 return (options == 0) ? setflag(status_flags::expanded) : *this;
 580         } else
 581                 return thisexpairseq(vp, overall_coeff);
 582 }
 583
 584 //////////
 585 // new virtual functions which can be overridden by derived classes
 586 //////////
 587
 588 // protected
 589
 590 /** Create an object of this type.
 591  *  This method works similar to a constructor.  It is useful because expairseq
 592  *  has (at least) two possible different semantics but we want to inherit
 593  *  methods thus avoiding code duplication.  Sometimes a method in expairseq
 594  *  has to create a new one of the same semantics, which cannot be done by a
 595  *  ctor because the name (add, mul,...) is unknown on the expaiseq level.  In
 596  *  order for this trick to work a derived class must of course override this
 597  *  definition. */
 598 ex expairseq::thisexpairseq(const epvector &v, const ex &oc) const
 599 {
 600         return expairseq(v,oc);
 601 }
 602
 603 ex expairseq::thisexpairseq(epvector *vp, const ex &oc) const
 604 {
 605         return expairseq(vp,oc);
 606 }
 607
 608 void expairseq::printpair(const print_context & c, const expair & p, unsigned upper_precedence) const
 609 {
 610         c.s << "[[";
 611         p.rest.print(c, precedence());
 612         c.s << ",";
 613         p.coeff.print(c, precedence());
 614         c.s << "]]";
 615 }
 616
 617 void expairseq::printseq(const print_context & c, char delim,
 618                          unsigned this_precedence,
 619                          unsigned upper_precedence) const
 620 {
 621         if (this_precedence <= upper_precedence)
 622                 c.s << "(";
 623         epvector::const_iterator it, it_last = seq.end() - 1;
 624         for (it=seq.begin(); it!=it_last; ++it) {
 625                 printpair(c, *it, this_precedence);
 626                 c.s << delim;
 627         }
 628         printpair(c, *it, this_precedence);
 629         if (!overall_coeff.is_equal(default_overall_coeff())) {
 630                 c.s << delim;
 631                 overall_coeff.print(c, this_precedence);
 632         }
 633
 634         if (this_precedence <= upper_precedence)
 635                 c.s << ")";
 636 }
 637
 638
 639 /** Form an expair from an ex, using the corresponding semantics.
 640  *  @see expairseq::recombine_pair_to_ex() */
 641 expair expairseq::split_ex_to_pair(const ex &e) const
 642 {
 643         return expair(e,_ex1);
 644 }
 645
 646
 647 expair expairseq::combine_ex_with_coeff_to_pair(const ex &e,
 648                                                 const ex &c) const
 649 {
 650         GINAC_ASSERT(is_exactly_a<numeric>(c));
 651
 652         return expair(e,c);
 653 }
 654
 655
 656 expair expairseq::combine_pair_with_coeff_to_pair(const expair &p,
 657                                                   const ex &c) const
 658 {
 659         GINAC_ASSERT(is_exactly_a<numeric>(p.coeff));
 660         GINAC_ASSERT(is_exactly_a<numeric>(c));
 661
 662         return expair(p.rest,ex_to<numeric>(p.coeff).mul_dyn(ex_to<numeric>(c)));
 663 }
 664
 665
 666 /** Form an ex out of an expair, using the corresponding semantics.
 667  *  @see expairseq::split_ex_to_pair() */
 668 ex expairseq::recombine_pair_to_ex(const expair &p) const
 669 {
 670         return lst(p.rest,p.coeff);
 671 }
 672
 673 bool expairseq::expair_needs_further_processing(epp it)
 674 {
 675 #if EXPAIRSEQ_USE_HASHTAB
 676         //#  error "FIXME: expair_needs_further_processing not yet implemented for hashtabs, sorry. A.F."
 677 #endif // EXPAIRSEQ_USE_HASHTAB
 678         return false;
 679 }
 680
 681 ex expairseq::default_overall_coeff(void) const
 682 {
 683         return _ex0;
 684 }
 685
 686 void expairseq::combine_overall_coeff(const ex &c)
 687 {
 688         GINAC_ASSERT(is_exactly_a<numeric>(overall_coeff));
 689         GINAC_ASSERT(is_exactly_a<numeric>(c));
 690         overall_coeff = ex_to<numeric>(overall_coeff).add_dyn(ex_to<numeric>(c));
 691 }
 692
 693 void expairseq::combine_overall_coeff(const ex &c1, const ex &c2)
 694 {
 695         GINAC_ASSERT(is_exactly_a<numeric>(overall_coeff));
 696         GINAC_ASSERT(is_exactly_a<numeric>(c1));
 697         GINAC_ASSERT(is_exactly_a<numeric>(c2));
 698         overall_coeff = ex_to<numeric>(overall_coeff).
 699                         add_dyn(ex_to<numeric>(c1).mul(ex_to<numeric>(c2)));
 700 }
 701
 702 bool expairseq::can_make_flat(const expair &p) const
 703 {
 704         return true;
 705 }
 706
 707
 708 //////////
 709 // non-virtual functions in this class
 710 //////////
 711
 712 void expairseq::construct_from_2_ex_via_exvector(const ex &lh, const ex &rh)
 713 {
 714         exvector v;
 715         v.reserve(2);
 716         v.push_back(lh);
 717         v.push_back(rh);
 718         construct_from_exvector(v);
 719 #if EXPAIRSEQ_USE_HASHTAB
 720         GINAC_ASSERT((hashtabsize==0)||(hashtabsize>=minhashtabsize));
 721         GINAC_ASSERT(hashtabsize==calc_hashtabsize(seq.size()));
 722 #endif // EXPAIRSEQ_USE_HASHTAB
 723 }
 724
 725 void expairseq::construct_from_2_ex(const ex &lh, const ex &rh)
 726 {
 727         if (ex_to<basic>(lh).tinfo()==this->tinfo()) {
 728                 if (ex_to<basic>(rh).tinfo()==this->tinfo()) {
 729 #if EXPAIRSEQ_USE_HASHTAB
 730                         unsigned totalsize = ex_to<expairseq>(lh).seq.size() +
 731                                              ex_to<expairseq>(rh).seq.size();
 732                         if (calc_hashtabsize(totalsize)!=0) {
 733                                 construct_from_2_ex_via_exvector(lh,rh);
 734                         } else {
 735 #endif // EXPAIRSEQ_USE_HASHTAB
 736                                 construct_from_2_expairseq(ex_to<expairseq>(lh),
 737                                                            ex_to<expairseq>(rh));
 738 #if EXPAIRSEQ_USE_HASHTAB
 739                         }
 740 #endif // EXPAIRSEQ_USE_HASHTAB
 741                         return;
 742                 } else {
 743 #if EXPAIRSEQ_USE_HASHTAB
 744                         unsigned totalsize = ex_to<expairseq>(lh).seq.size()+1;
 745                         if (calc_hashtabsize(totalsize)!=0) {
 746                                 construct_from_2_ex_via_exvector(lh, rh);
 747                         } else {
 748 #endif // EXPAIRSEQ_USE_HASHTAB
 749                                 construct_from_expairseq_ex(ex_to<expairseq>(lh), rh);
 750 #if EXPAIRSEQ_USE_HASHTAB
 751                         }
 752 #endif // EXPAIRSEQ_USE_HASHTAB
 753                         return;
 754                 }
 755         } else if (ex_to<basic>(rh).tinfo()==this->tinfo()) {
 756 #if EXPAIRSEQ_USE_HASHTAB
 757                 unsigned totalsize=ex_to<expairseq>(rh).seq.size()+1;
 758                 if (calc_hashtabsize(totalsize)!=0) {
 759                         construct_from_2_ex_via_exvector(lh,rh);
 760                 } else {
 761 #endif // EXPAIRSEQ_USE_HASHTAB
 762                         construct_from_expairseq_ex(ex_to<expairseq>(rh),lh);
 763 #if EXPAIRSEQ_USE_HASHTAB
 764                 }
 765 #endif // EXPAIRSEQ_USE_HASHTAB
 766                 return;
 767         }
 768
 769 #if EXPAIRSEQ_USE_HASHTAB
 770         if (calc_hashtabsize(2)!=0) {
 771                 construct_from_2_ex_via_exvector(lh,rh);
 772                 return;
 773         }
 774         hashtabsize = 0;
 775 #endif // EXPAIRSEQ_USE_HASHTAB
 776
 777         if (is_ex_exactly_of_type(lh,numeric)) {
 778                 if (is_ex_exactly_of_type(rh,numeric)) {
 779                         combine_overall_coeff(lh);
 780                         combine_overall_coeff(rh);
 781                 } else {
 782                         combine_overall_coeff(lh);
 783                         seq.push_back(split_ex_to_pair(rh));
 784                 }
 785         } else {
 786                 if (is_ex_exactly_of_type(rh,numeric)) {
 787                         combine_overall_coeff(rh);
 788                         seq.push_back(split_ex_to_pair(lh));
 789                 } else {
 790                         expair p1 = split_ex_to_pair(lh);
 791                         expair p2 = split_ex_to_pair(rh);
 792
 793                         int cmpval = p1.rest.compare(p2.rest);
 794                         if (cmpval==0) {
 795                                 p1.coeff = ex_to<numeric>(p1.coeff).add_dyn(ex_to<numeric>(p2.coeff));
 796                                 if (!ex_to<numeric>(p1.coeff).is_zero()) {
 797                                         // no further processing is necessary, since this
 798                                         // one element will usually be recombined in eval()
 799                                         seq.push_back(p1);
 800                                 }
 801                         } else {
 802                                 seq.reserve(2);
 803                                 if (cmpval<0) {
 804                                         seq.push_back(p1);
 805                                         seq.push_back(p2);
 806                                 } else {
 807                                         seq.push_back(p2);
 808                                         seq.push_back(p1);
 809                                 }
 810                         }
 811                 }
 812         }
 813 }
 814
 815 void expairseq::construct_from_2_expairseq(const expairseq &s1,
 816                                                                                    const expairseq &s2)
 817 {
 818         combine_overall_coeff(s1.overall_coeff);
 819         combine_overall_coeff(s2.overall_coeff);
 820
 821         epvector::const_iterator first1 = s1.seq.begin();
 822         epvector::const_iterator last1 = s1.seq.end();
 823         epvector::const_iterator first2 = s2.seq.begin();
 824         epvector::const_iterator last2 = s2.seq.end();
 825
 826         seq.reserve(s1.seq.size()+s2.seq.size());
 827
 828         bool needs_further_processing=false;
 829
 830         while (first1!=last1 && first2!=last2) {
 831                 int cmpval = (*first1).rest.compare((*first2).rest);
 832                 if (cmpval==0) {
 833                         // combine terms
 834                         const numeric &newcoeff = ex_to<numeric>(first1->coeff).
 835                                                    add(ex_to<numeric>(first2->coeff));
 836                         if (!newcoeff.is_zero()) {
 837                                 seq.push_back(expair(first1->rest,newcoeff));
 838                                 if (expair_needs_further_processing(seq.end()-1)) {
 839                                         needs_further_processing = true;
 840                                 }
 841                         }
 842                         ++first1;
 843                         ++first2;
 844                 } else if (cmpval<0) {
 845                         seq.push_back(*first1);
 846                         ++first1;
 847                 } else {
 848                         seq.push_back(*first2);
 849                         ++first2;
 850                 }
 851         }
 852
 853         while (first1!=last1) {
 854                 seq.push_back(*first1);
 855                 ++first1;
 856         }
 857         while (first2!=last2) {
 858                 seq.push_back(*first2);
 859                 ++first2;
 860         }
 861
 862         if (needs_further_processing) {
 863                 epvector v = seq;
 864                 seq.clear();
 865                 construct_from_epvector(v);
 866         }
 867 }
 868
 869 void expairseq::construct_from_expairseq_ex(const expairseq &s,
 870                                                                                         const ex &e)
 871 {
 872         combine_overall_coeff(s.overall_coeff);
 873         if (is_ex_exactly_of_type(e,numeric)) {
 874                 combine_overall_coeff(e);
 875                 seq = s.seq;
 876                 return;
 877         }
 878
 879         epvector::const_iterator first = s.seq.begin();
 880         epvector::const_iterator last = s.seq.end();
 881         expair p = split_ex_to_pair(e);
 882
 883         seq.reserve(s.seq.size()+1);
 884         bool p_pushed = false;
 885
 886         bool needs_further_processing=false;
 887
 888         // merge p into s.seq
 889         while (first!=last) {
 890                 int cmpval = (*first).rest.compare(p.rest);
 891                 if (cmpval==0) {
 892                         // combine terms
 893                         const numeric &newcoeff = ex_to<numeric>(first->coeff).
 894                                                    add(ex_to<numeric>(p.coeff));
 895                         if (!newcoeff.is_zero()) {
 896                                 seq.push_back(expair(first->rest,newcoeff));
 897                                 if (expair_needs_further_processing(seq.end()-1))
 898                                         needs_further_processing = true;
 899                         }
 900                         ++first;
 901                         p_pushed = true;
 902                         break;
 903                 } else if (cmpval<0) {
 904                         seq.push_back(*first);
 905                         ++first;
 906                 } else {
 907                         seq.push_back(p);
 908                         p_pushed = true;
 909                         break;
 910                 }
 911         }
 912
 913         if (p_pushed) {
 914                 // while loop exited because p was pushed, now push rest of s.seq
 915                 while (first!=last) {
 916                         seq.push_back(*first);
 917                         ++first;
 918                 }
 919         } else {
 920                 // while loop exited because s.seq was pushed, now push p
 921                 seq.push_back(p);
 922         }
 923
 924         if (needs_further_processing) {
 925                 epvector v = seq;
 926                 seq.clear();
 927                 construct_from_epvector(v);
 928         }
 929 }
 930
 931 void expairseq::construct_from_exvector(const exvector &v)
 932 {
 933         // simplifications: +(a,+(b,c),d) -> +(a,b,c,d) (associativity)
 934         //                  +(d,b,c,a) -> +(a,b,c,d) (canonicalization)
 935         //                  +(...,x,*(x,c1),*(x,c2)) -> +(...,*(x,1+c1+c2)) (c1, c2 numeric())
 936         //                  (same for (+,*) -> (*,^)
 937
 938         make_flat(v);
 939 #if EXPAIRSEQ_USE_HASHTAB
 940         combine_same_terms();
 941 #else
 942         canonicalize();
 943         combine_same_terms_sorted_seq();
 944 #endif // EXPAIRSEQ_USE_HASHTAB
 945 }
 946
 947 void expairseq::construct_from_epvector(const epvector &v)
 948 {
 949         // simplifications: +(a,+(b,c),d) -> +(a,b,c,d) (associativity)
 950         //                  +(d,b,c,a) -> +(a,b,c,d) (canonicalization)
 951         //                  +(...,x,*(x,c1),*(x,c2)) -> +(...,*(x,1+c1+c2)) (c1, c2 numeric())
 952         //                  (same for (+,*) -> (*,^)
 953
 954         make_flat(v);
 955 #if EXPAIRSEQ_USE_HASHTAB
 956         combine_same_terms();
 957 #else
 958         canonicalize();
 959         combine_same_terms_sorted_seq();
 960 #endif // EXPAIRSEQ_USE_HASHTAB
 961 }
 962
 963 /** Combine this expairseq with argument exvector.
 964  *  It cares for associativity as well as for special handling of numerics. */
 965 void expairseq::make_flat(const exvector &v)
 966 {
 967         exvector::const_iterator cit;
 968
 969         // count number of operands which are of same expairseq derived type
 970         // and their cumulative number of operands
 971         int nexpairseqs = 0;
 972         int noperands = 0;
 973
 974         cit = v.begin();
 975         while (cit!=v.end()) {
 976                 if (ex_to<basic>(*cit).tinfo()==this->tinfo()) {
 977                         ++nexpairseqs;
 978                         noperands += ex_to<expairseq>(*cit).seq.size();
 979                 }
 980                 ++cit;
 981         }
 982
 983         // reserve seq and coeffseq which will hold all operands
 984         seq.reserve(v.size()+noperands-nexpairseqs);
 985
 986         // copy elements and split off numerical part
 987         cit = v.begin();
 988         while (cit!=v.end()) {
 989                 if (ex_to<basic>(*cit).tinfo()==this->tinfo()) {
 990                         const expairseq &subseqref = ex_to<expairseq>(*cit);
 991                         combine_overall_coeff(subseqref.overall_coeff);
 992                         epvector::const_iterator cit_s = subseqref.seq.begin();
 993                         while (cit_s!=subseqref.seq.end()) {
 994                                 seq.push_back(*cit_s);
 995                                 ++cit_s;
 996                         }
 997                 } else {
 998                         if (is_ex_exactly_of_type(*cit,numeric))
 999                                 combine_overall_coeff(*cit);
1000                         else
1001                                 seq.push_back(split_ex_to_pair(*cit));
1002                 }
1003                 ++cit;
1004         }
1005 }
1006
1007 /** Combine this expairseq with argument epvector.
1008  *  It cares for associativity as well as for special handling of numerics. */
1009 void expairseq::make_flat(const epvector &v)
1010 {
1011         epvector::const_iterator cit;
1012
1013         // count number of operands which are of same expairseq derived type
1014         // and their cumulative number of operands
1015         int nexpairseqs = 0;
1016         int noperands = 0;
1017
1018         cit = v.begin();
1019         while (cit!=v.end()) {
1020                 if (ex_to<basic>(cit->rest).tinfo()==this->tinfo()) {
1021                         ++nexpairseqs;
1022                         noperands += ex_to<expairseq>(cit->rest).seq.size();
1023                 }
1024                 ++cit;
1025         }
1026
1027         // reserve seq and coeffseq which will hold all operands
1028         seq.reserve(v.size()+noperands-nexpairseqs);
1029
1030         // copy elements and split off numerical part
1031         cit = v.begin();
1032         while (cit!=v.end()) {
1033                 if (ex_to<basic>(cit->rest).tinfo()==this->tinfo() &&
1034                     this->can_make_flat(*cit)) {
1035                         const expairseq &subseqref = ex_to<expairseq>(cit->rest);
1036                         combine_overall_coeff(ex_to<numeric>(subseqref.overall_coeff),
1037                                                             ex_to<numeric>(cit->coeff));
1038                         epvector::const_iterator cit_s = subseqref.seq.begin();
1039                         while (cit_s!=subseqref.seq.end()) {
1040                                 seq.push_back(expair(cit_s->rest,
1041                                                      ex_to<numeric>(cit_s->coeff).mul_dyn(ex_to<numeric>(cit->coeff))));
1042                                 //seq.push_back(combine_pair_with_coeff_to_pair(*cit_s,
1043                                 //                                              (*cit).coeff));
1044                                 ++cit_s;
1045                         }
1046                 } else {
1047                         if (cit->is_canonical_numeric())
1048                                 combine_overall_coeff(cit->rest);
1049                         else
1050                                 seq.push_back(*cit);
1051                 }
1052                 ++cit;
1053         }
1054 }
1055
1056 /** Brings this expairseq into a sorted (canonical) form. */
1057 void expairseq::canonicalize(void)
1058 {
1059         sort(seq.begin(), seq.end(), expair_is_less());
1060 }
1061
1062
1063 /** Compact a presorted expairseq by combining all matching expairs to one
1064  *  each.  On an add object, this is responsible for 2*x+3*x+y -> 5*x+y, for
1065  *  instance. */
1066 void expairseq::combine_same_terms_sorted_seq(void)
1067 {
1068         bool needs_further_processing = false;
1069
1070         if (seq.size()>1) {
1071                 epvector::iterator itin1 = seq.begin();
1072                 epvector::iterator itin2 = itin1+1;
1073                 epvector::iterator itout = itin1;
1074                 epvector::iterator last = seq.end();
1075                 // must_copy will be set to true the first time some combination is
1076                 // possible from then on the sequence has changed and must be compacted
1077                 bool must_copy = false;
1078                 while (itin2!=last) {
1079                         if (itin1->rest.compare(itin2->rest)==0) {
1080                                 itin1->coeff = ex_to<numeric>(itin1->coeff).
1081                                                add_dyn(ex_to<numeric>(itin2->coeff));
1082                                 if (expair_needs_further_processing(itin1))
1083                                         needs_further_processing = true;
1084                                 must_copy = true;
1085                         } else {
1086                                 if (!ex_to<numeric>(itin1->coeff).is_zero()) {
1087                                         if (must_copy)
1088                                                 *itout = *itin1;
1089                                         ++itout;
1090                                 }
1091                                 itin1 = itin2;
1092                         }
1093                         ++itin2;
1094                 }
1095                 if (!ex_to<numeric>(itin1->coeff).is_zero()) {
1096                         if (must_copy)
1097                                 *itout = *itin1;
1098                         ++itout;
1099                 }
1100                 if (itout!=last)
1101                         seq.erase(itout,last);
1102         }
1103
1104         if (needs_further_processing) {
1105                 epvector v = seq;
1106                 seq.clear();
1107                 construct_from_epvector(v);
1108         }
1109 }
1110
1111 #if EXPAIRSEQ_USE_HASHTAB
1112
1113 unsigned expairseq::calc_hashtabsize(unsigned sz) const
1114 {
1115         unsigned size;
1116         unsigned nearest_power_of_2 = 1 << log2(sz);
1117         // if (nearest_power_of_2 < maxhashtabsize/hashtabfactor) {
1118         //  size = nearest_power_of_2*hashtabfactor;
1119         size = nearest_power_of_2/hashtabfactor;
1120         if (size<minhashtabsize)
1121                 return 0;
1122         GINAC_ASSERT(hashtabsize<=0x8000000U); // really max size due to 31 bit hashing
1123         // hashtabsize must be a power of 2
1124         GINAC_ASSERT((1U << log2(size))==size);
1125         return size;
1126 }
1127
1128 unsigned expairseq::calc_hashindex(const ex &e) const
1129 {
1130         // calculate hashindex
1131         unsigned hash = e.gethash();
1132         unsigned hashindex;
1133         if (is_a_numeric_hash(hash)) {
1134                 hashindex = hashmask;
1135         } else {
1136                 hashindex = hash &hashmask;
1137                 // last hashtab entry is reserved for numerics
1138                 if (hashindex==hashmask) hashindex = 0;
1139         }
1140         GINAC_ASSERT(hashindex>=0);
1141         GINAC_ASSERT((hashindex<hashtabsize)||(hashtabsize==0));
1142         return hashindex;
1143 }
1144
1145 void expairseq::shrink_hashtab(void)
1146 {
1147         unsigned new_hashtabsize;
1148         while (hashtabsize!=(new_hashtabsize=calc_hashtabsize(seq.size()))) {
1149                 GINAC_ASSERT(new_hashtabsize<hashtabsize);
1150                 if (new_hashtabsize==0) {
1151                         hashtab.clear();
1152                         hashtabsize = 0;
1153                         canonicalize();
1154                         return;
1155                 }
1156
1157                 // shrink by a factor of 2
1158                 unsigned half_hashtabsize = hashtabsize/2;
1159                 for (unsigned i=0; i<half_hashtabsize-1; ++i)
1160                         hashtab[i].merge(hashtab[i+half_hashtabsize],epp_is_less());
1161                 // special treatment for numeric hashes
1162                 hashtab[0].merge(hashtab[half_hashtabsize-1],epp_is_less());
1163                 hashtab[half_hashtabsize-1] = hashtab[hashtabsize-1];
1164                 hashtab.resize(half_hashtabsize);
1165                 hashtabsize = half_hashtabsize;
1166                 hashmask = hashtabsize-1;
1167         }
1168 }
1169
1170 void expairseq::remove_hashtab_entry(epvector::const_iterator element)
1171 {
1172         if (hashtabsize==0)
1173                 return; // nothing to do
1174
1175         // calculate hashindex of element to be deleted
1176         unsigned hashindex = calc_hashindex((*element).rest);
1177
1178         // find it in hashtab and remove it
1179         epplist &eppl = hashtab[hashindex];
1180         epplist::iterator epplit = eppl.begin();
1181         bool erased = false;
1182         while (epplit!=eppl.end()) {
1183                 if (*epplit == element) {
1184                         eppl.erase(epplit);
1185                         erased = true;
1186                         break;
1187                 }
1188                 ++epplit;
1189         }
1190         if (!erased) {
1191                 std::cout << "tried to erase " << element-seq.begin() << std::endl;
1192                 std::cout << "size " << seq.end()-seq.begin() << std::endl;
1193
1194                 unsigned hashindex = calc_hashindex(element->rest);
1195                 epplist &eppl = hashtab[hashindex];
1196                 epplist::iterator epplit = eppl.begin();
1197                 bool erased = false;
1198                 while (epplit!=eppl.end()) {
1199                         if (*epplit == element) {
1200                                 eppl.erase(epplit);
1201                                 erased = true;
1202                                 break;
1203                         }
1204                         ++epplit;
1205                 }
1206                 GINAC_ASSERT(erased);
1207         }
1208         GINAC_ASSERT(erased);
1209 }
1210
1211 void expairseq::move_hashtab_entry(epvector::const_iterator oldpos,
1212                                    epvector::iterator newpos)
1213 {
1214         GINAC_ASSERT(hashtabsize!=0);
1215
1216         // calculate hashindex of element which was moved
1217         unsigned hashindex=calc_hashindex((*newpos).rest);
1218
1219         // find it in hashtab and modify it
1220         epplist &eppl = hashtab[hashindex];
1221         epplist::iterator epplit = eppl.begin();
1222         while (epplit!=eppl.end()) {
1223                 if (*epplit == oldpos) {
1224                         *epplit = newpos;
1225                         break;
1226                 }
1227                 ++epplit;
1228         }
1229         GINAC_ASSERT(epplit!=eppl.end());
1230 }
1231
1232 void expairseq::sorted_insert(epplist &eppl, epvector::const_iterator elem)
1233 {
1234         epplist::const_iterator current = eppl.begin();
1235         while ((current!=eppl.end()) && ((*current)->is_less(*elem))) {
1236                 ++current;
1237         }
1238         eppl.insert(current,elem);
1239 }
1240
1241 void expairseq::build_hashtab_and_combine(epvector::iterator &first_numeric,
1242                                           epvector::iterator &last_non_zero,
1243                                           std::vector<bool> &touched,
1244                                           unsigned &number_of_zeroes)
1245 {
1246         epp current = seq.begin();
1247
1248         while (current!=first_numeric) {
1249                 if (is_ex_exactly_of_type(current->rest,numeric)) {
1250                         --first_numeric;
1251                         iter_swap(current,first_numeric);
1252                 } else {
1253                         // calculate hashindex
1254                         unsigned currenthashindex = calc_hashindex(current->rest);
1255
1256                         // test if there is already a matching expair in the hashtab-list
1257                         epplist &eppl=hashtab[currenthashindex];
1258                         epplist::iterator epplit = eppl.begin();
1259                         while (epplit!=eppl.end()) {
1260                                 if (current->rest.is_equal((*epplit)->rest))
1261                                         break;
1262                                 ++epplit;
1263                         }
1264                         if (epplit==eppl.end()) {
1265                                 // no matching expair found, append this to end of list
1266                                 sorted_insert(eppl,current);
1267                                 ++current;
1268                         } else {
1269                                 // epplit points to a matching expair, combine it with current
1270                                 (*epplit)->coeff = ex_to<numeric>((*epplit)->coeff).
1271                                                    add_dyn(ex_to<numeric>(current->coeff));
1272
1273                                 // move obsolete current expair to end by swapping with last_non_zero element
1274                                 // if this was a numeric, it is swapped with the expair before first_numeric
1275                                 iter_swap(current,last_non_zero);
1276                                 --first_numeric;
1277                                 if (first_numeric!=last_non_zero) iter_swap(first_numeric,current);
1278                                 --last_non_zero;
1279                                 ++number_of_zeroes;
1280                                 // test if combined term has coeff 0 and can be removed is done later
1281                                 touched[(*epplit)-seq.begin()] = true;
1282                         }
1283                 }
1284         }
1285 }
1286
1287 void expairseq::drop_coeff_0_terms(epvector::iterator &first_numeric,
1288                                    epvector::iterator &last_non_zero,
1289                                    std::vector<bool> &touched,
1290                                    unsigned &number_of_zeroes)
1291 {
1292         // move terms with coeff 0 to end and remove them from hashtab
1293         // check only those elements which have been touched
1294         epp current = seq.begin();
1295         unsigned i = 0;
1296         while (current!=first_numeric) {
1297                 if (!touched[i]) {
1298                         ++current;
1299                         ++i;
1300                 } else if (!ex_to<numeric>((*current).coeff).is_zero()) {
1301                         ++current;
1302                         ++i;
1303                 } else {
1304                         remove_hashtab_entry(current);
1305
1306                         // move element to the end, unless it is already at the end
1307                         if (current!=last_non_zero) {
1308                                 iter_swap(current,last_non_zero);
1309                                 --first_numeric;
1310                                 bool numeric_swapped = first_numeric!=last_non_zero;
1311                                 if (numeric_swapped)
1312                                         iter_swap(first_numeric,current);
1313                                 epvector::iterator changed_entry;
1314
1315                                 if (numeric_swapped)
1316                                         changed_entry = first_numeric;
1317                                 else
1318                                         changed_entry = last_non_zero;
1319
1320                                 --last_non_zero;
1321                                 ++number_of_zeroes;
1322
1323                                 if (first_numeric!=current) {
1324
1325                                         // change entry in hashtab which referred to first_numeric or last_non_zero to current
1326                                         move_hashtab_entry(changed_entry,current);
1327                                         touched[current-seq.begin()] = touched[changed_entry-seq.begin()];
1328                                 }
1329                         } else {
1330                                 --first_numeric;
1331                                 --last_non_zero;
1332                                 ++number_of_zeroes;
1333                         }
1334                 }
1335         }
1336         GINAC_ASSERT(i==current-seq.begin());
1337 }
1338
1339 /** True if one of the coeffs vanishes, otherwise false.
1340  *  This would be an invariant violation, so this should only be used for
1341  *  debugging purposes. */
1342 bool expairseq::has_coeff_0(void) const
1343 {
1344         epvector::const_iterator i = seq.begin(), end = seq.end();
1345         while (i != end) {
1346                 if (i->coeff.is_zero())
1347                         return true;
1348                 ++i;
1349         }
1350         return false;
1351 }
1352
1353 void expairseq::add_numerics_to_hashtab(epvector::iterator first_numeric,
1354                                                                                 epvector::const_iterator last_non_zero)
1355 {
1356         if (first_numeric == seq.end()) return; // no numerics
1357
1358         epvector::const_iterator current = first_numeric, last = last_non_zero + 1;
1359         while (current != last) {
1360                 sorted_insert(hashtab[hashmask], current);
1361                 ++current;
1362         }
1363 }
1364
1365 void expairseq::combine_same_terms(void)
1366 {
1367         // combine same terms, drop term with coeff 0, move numerics to end
1368
1369         // calculate size of hashtab
1370         hashtabsize = calc_hashtabsize(seq.size());
1371
1372         // hashtabsize is a power of 2
1373         hashmask = hashtabsize-1;
1374
1375         // allocate hashtab
1376         hashtab.clear();
1377         hashtab.resize(hashtabsize);
1378
1379         if (hashtabsize==0) {
1380                 canonicalize();
1381                 combine_same_terms_sorted_seq();
1382                 GINAC_ASSERT(!has_coeff_0());
1383                 return;
1384         }
1385
1386         // iterate through seq, move numerics to end,
1387         // fill hashtab and combine same terms
1388         epvector::iterator first_numeric = seq.end();
1389         epvector::iterator last_non_zero = seq.end()-1;
1390
1391         unsigned num = seq.size();
1392         std::vector<bool> touched(num);
1393
1394         unsigned number_of_zeroes = 0;
1395
1396         GINAC_ASSERT(!has_coeff_0());
1397         build_hashtab_and_combine(first_numeric,last_non_zero,touched,number_of_zeroes);
1398
1399         // there should not be any terms with coeff 0 from the beginning,
1400         // so it should be safe to skip this step
1401         if (number_of_zeroes!=0) {
1402                 drop_coeff_0_terms(first_numeric,last_non_zero,touched,number_of_zeroes);
1403         }
1404
1405         add_numerics_to_hashtab(first_numeric,last_non_zero);
1406
1407         // pop zero elements
1408         for (unsigned i=0; i<number_of_zeroes; ++i) {
1409                 seq.pop_back();
1410         }
1411
1412         // shrink hashtabsize to calculated value
1413         GINAC_ASSERT(!has_coeff_0());
1414
1415         shrink_hashtab();
1416
1417         GINAC_ASSERT(!has_coeff_0());
1418 }
1419
1420 #endif // EXPAIRSEQ_USE_HASHTAB
1421
1422 /** Check if this expairseq is in sorted (canonical) form.  Useful mainly for
1423  *  debugging or in assertions since being sorted is an invariance. */
1424 bool expairseq::is_canonical() const
1425 {
1426         if (seq.size() <= 1)
1427                 return 1;
1428
1429 #if EXPAIRSEQ_USE_HASHTAB
1430         if (hashtabsize > 0) return 1; // not canoncalized
1431 #endif // EXPAIRSEQ_USE_HASHTAB
1432
1433         epvector::const_iterator it = seq.begin(), itend = seq.end();
1434         epvector::const_iterator it_last = it;
1435         for (++it; it!=itend; it_last=it, ++it) {
1436                 if (!(it_last->is_less(*it) || it_last->is_equal(*it))) {
1437                         if (!is_ex_exactly_of_type(it_last->rest,numeric) ||
1438                                 !is_ex_exactly_of_type(it->rest,numeric)) {
1439                                 // double test makes it easier to set a breakpoint...
1440                                 if (!is_ex_exactly_of_type(it_last->rest,numeric) ||
1441                                         !is_ex_exactly_of_type(it->rest,numeric)) {
1442                                         printpair(std::clog, *it_last, 0);
1443                                         std::clog << ">";
1444                                         printpair(std::clog, *it, 0);
1445                                         std::clog << "\n";
1446                                         std::clog << "pair1:" << std::endl;
1447                                         it_last->rest.print(print_tree(std::clog));
1448                                         it_last->coeff.print(print_tree(std::clog));
1449                                         std::clog << "pair2:" << std::endl;
1450                                         it->rest.print(print_tree(std::clog));
1451                                         it->coeff.print(print_tree(std::clog));
1452                                         return 0;
1453                                 }
1454                         }
1455                 }
1456         }
1457         return 1;
1458 }
1459
1460
1461 /** Member-wise expand the expairs in this sequence.
1462  *
1463  *  @see expairseq::expand()
1464  *  @return pointer to epvector containing expanded pairs or zero pointer,
1465  *  if no members were changed. */
1466 epvector * expairseq::expandchildren(unsigned options) const
1467 {
1468         epvector::const_iterator last = seq.end();
1469         epvector::const_iterator cit = seq.begin();
1470         while (cit!=last) {
1471                 const ex &expanded_ex = cit->rest.expand(options);
1472                 if (!are_ex_trivially_equal(cit->rest,expanded_ex)) {
1473
1474                         // something changed, copy seq, eval and return it
1475                         epvector *s = new epvector;
1476                         s->reserve(seq.size());
1477
1478                         // copy parts of seq which are known not to have changed
1479                         epvector::const_iterator cit2 = seq.begin();
1480                         while (cit2!=cit) {
1481                                 s->push_back(*cit2);
1482                                 ++cit2;
1483                         }
1484                         // copy first changed element
1485                         s->push_back(combine_ex_with_coeff_to_pair(expanded_ex,
1486                                                                    cit2->coeff));
1487                         ++cit2;
1488                         // copy rest
1489                         while (cit2!=last) {
1490                                 s->push_back(combine_ex_with_coeff_to_pair(cit2->rest.expand(options),
1491                                                                            cit2->coeff));
1492                                 ++cit2;
1493                         }
1494                         return s;
1495                 }
1496                 ++cit;
1497         }
1498
1499         return 0; // signalling nothing has changed
1500 }
1501
1502
1503 /** Member-wise evaluate the expairs in this sequence.
1504  *
1505  *  @see expairseq::eval()
1506  *  @return pointer to epvector containing evaluated pairs or zero pointer,
1507  *  if no members were changed. */
1508 epvector * expairseq::evalchildren(int level) const
1509 {
1510         // returns a NULL pointer if nothing had to be evaluated
1511         // returns a pointer to a newly created epvector otherwise
1512         // (which has to be deleted somewhere else)
1513
1514         if (level==1)
1515                 return 0;
1516
1517         if (level == -max_recursion_level)
1518                 throw(std::runtime_error("max recursion level reached"));
1519
1520         --level;
1521         epvector::const_iterator last = seq.end();
1522         epvector::const_iterator cit = seq.begin();
1523         while (cit!=last) {
1524                 const ex &evaled_ex = cit->rest.eval(level);
1525                 if (!are_ex_trivially_equal(cit->rest,evaled_ex)) {
1526
1527                         // something changed, copy seq, eval and return it
1528                         epvector *s = new epvector;
1529                         s->reserve(seq.size());
1530
1531                         // copy parts of seq which are known not to have changed
1532                         epvector::const_iterator cit2=seq.begin();
1533                         while (cit2!=cit) {
1534                                 s->push_back(*cit2);
1535                                 ++cit2;
1536                         }
1537                         // copy first changed element
1538                         s->push_back(combine_ex_with_coeff_to_pair(evaled_ex,
1539                                                                    cit2->coeff));
1540                         ++cit2;
1541                         // copy rest
1542                         while (cit2!=last) {
1543                                 s->push_back(combine_ex_with_coeff_to_pair(cit2->rest.eval(level),
1544                                                                            cit2->coeff));
1545                                 ++cit2;
1546                         }
1547                         return s;
1548                 }
1549                 ++cit;
1550         }
1551
1552         return 0; // signalling nothing has changed
1553 }
1554
1555
1556 /** Member-wise substitute in this sequence.
1557  *
1558  *  @see expairseq::subs()
1559  *  @return pointer to epvector containing pairs after application of subs,
1560  *    or NULL pointer if no members were changed. */
1561 epvector * expairseq::subschildren(const lst &ls, const lst &lr, bool no_pattern) const
1562 {
1563         GINAC_ASSERT(ls.nops()==lr.nops());
1564
1565         // The substitution is "complex" when any of the objects to be substituted
1566         // is a product or power. In this case we have to recombine the pairs
1567         // because the numeric coefficients may be part of the search pattern.
1568         bool complex_subs = false;
1569         for (unsigned i=0; i<ls.nops(); ++i)
1570                 if (is_ex_exactly_of_type(ls.op(i), mul) || is_ex_exactly_of_type(ls.op(i), power)) {
1571                         complex_subs = true;
1572                         break;
1573                 }
1574
1575         if (complex_subs) {
1576
1577                 // Substitute in the recombined pairs
1578                 epvector::const_iterator cit = seq.begin(), last = seq.end();
1579                 while (cit != last) {
1580
1581                         const ex &orig_ex = recombine_pair_to_ex(*cit);
1582                         const ex &subsed_ex = orig_ex.subs(ls, lr, no_pattern);
1583                         if (!are_ex_trivially_equal(orig_ex, subsed_ex)) {
1584
1585                                 // Something changed, copy seq, subs and return it
1586                                 epvector *s = new epvector;
1587                                 s->reserve(seq.size());
1588
1589                                 // Copy parts of seq which are known not to have changed
1590                                 s->insert(s->begin(), seq.begin(), cit);
1591
1592                                 // Copy first changed element
1593                                 s->push_back(split_ex_to_pair(subsed_ex));
1594                                 ++cit;
1595
1596                                 // Copy rest
1597                                 while (cit != last) {
1598                                         s->push_back(split_ex_to_pair(recombine_pair_to_ex(*cit).subs(ls, lr, no_pattern)));
1599                                         ++cit;
1600                                 }
1601                                 return s;
1602                         }
1603
1604                         ++cit;
1605                 }
1606
1607         } else {
1608
1609                 // Substitute only in the "rest" part of the pairs
1610                 epvector::const_iterator cit = seq.begin(), last = seq.end();
1611                 while (cit != last) {
1612
1613                         const ex &subsed_ex = cit->rest.subs(ls, lr, no_pattern);
1614                         if (!are_ex_trivially_equal(cit->rest, subsed_ex)) {
1615
1616                                 // Something changed, copy seq, subs and return it
1617                                 epvector *s = new epvector;
1618                                 s->reserve(seq.size());
1619
1620                                 // Copy parts of seq which are known not to have changed
1621                                 s->insert(s->begin(), seq.begin(), cit);
1622
1623                                 // Copy first changed element
1624                                 s->push_back(combine_ex_with_coeff_to_pair(subsed_ex, cit->coeff));
1625                                 ++cit;
1626
1627                                 // Copy rest
1628                                 while (cit != last) {
1629                                         s->push_back(combine_ex_with_coeff_to_pair(cit->rest.subs(ls, lr, no_pattern),
1630                                                                                    cit->coeff));
1631                                         ++cit;
1632                                 }
1633                                 return s;
1634                         }
1635
1636                         ++cit;
1637                 }
1638         }
1639
1640         // Nothing has changed
1641         return NULL;
1642 }
1643
1644 //////////
1645 // static member variables
1646 //////////
1647
1648 #if EXPAIRSEQ_USE_HASHTAB
1649 unsigned expairseq::maxhashtabsize = 0x4000000U;
1650 unsigned expairseq::minhashtabsize = 0x1000U;
1651 unsigned expairseq::hashtabfactor = 1;
1652 #endif // EXPAIRSEQ_USE_HASHTAB
1653
1654 } // namespace GiNaC