ginac/expairseq.cpp

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