ginac/expairseq.cpp

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