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