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