ginac/expairseq.cpp

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