ginac/expairseq.cpp

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