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