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