ginac/expairseq.cpp

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