ginac/expairseq.cpp

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