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