1 /** @file expairseq.cpp
3 * Implementation of sequences of expression pairs. */
6 * GiNaC Copyright (C) 1999 Johannes Gutenberg University Mainz, Germany
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.
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.
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
28 #include "expairseq.h"
33 #ifndef NO_GINAC_NAMESPACE
35 #endif // ndef NO_GINAC_NAMESPACE
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
48 bool operator()(epp const & lh, epp const & rh) const
50 return (*lh).is_less(*rh);
55 // default constructor, destructor, copy constructor assignment operator and helpers
60 expairseq::expairseq(expairseq const & other)
62 debugmsg("expairseq copy constructor",LOGLEVEL_CONSTRUCT);
66 expairseq const & expairseq::operator=(expairseq const & other)
68 debugmsg("expairseq operator=",LOGLEVEL_ASSIGNMENT);
78 void expairseq::copy(expairseq const & other)
82 overall_coeff=other.overall_coeff;
83 #ifdef EXPAIRSEQ_USE_HASHTAB
85 hashtabsize=other.hashtabsize;
87 hashmask=other.hashmask;
88 hashtab.resize(hashtabsize);
89 epvector::const_iterator osb=other.seq.begin();
90 for (unsigned i=0; i<hashtabsize; ++i) {
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));
100 #endif // def EXPAIRSEQ_USE_HASHTAB
104 // other constructors
107 expairseq::expairseq(ex const & lh, ex const & rh) : basic(TINFO_expairseq)
109 debugmsg("expairseq constructor from ex,ex",LOGLEVEL_CONSTRUCT);
110 construct_from_2_ex(lh,rh);
111 GINAC_ASSERT(is_canonical());
114 expairseq::expairseq(exvector const & v) : basic(TINFO_expairseq)
116 debugmsg("expairseq constructor from exvector",LOGLEVEL_CONSTRUCT);
117 construct_from_exvector(v);
118 GINAC_ASSERT(is_canonical());
122 expairseq::expairseq(epvector const & v, bool do_not_canonicalize) :
123 basic(TINFO_expairseq)
125 debugmsg("expairseq constructor from epvector",LOGLEVEL_CONSTRUCT);
126 if (do_not_canonicalize) {
128 #ifdef EXPAIRSEQ_USE_HASHTAB
129 combine_same_terms(); // to build hashtab
130 #endif // def EXPAIRSEQ_USE_HASHTAB
132 construct_from_epvector(v);
134 GINAC_ASSERT(is_canonical());
138 expairseq::expairseq(epvector const & v, ex const & oc) :
139 basic(TINFO_expairseq), overall_coeff(oc)
141 debugmsg("expairseq constructor from epvector,ex",LOGLEVEL_CONSTRUCT);
142 construct_from_epvector(v);
143 GINAC_ASSERT(is_canonical());
146 expairseq::expairseq(epvector * vp, ex const & oc) :
147 basic(TINFO_expairseq), overall_coeff(oc)
149 debugmsg("expairseq constructor from epvector *,ex",LOGLEVEL_CONSTRUCT);
151 construct_from_epvector(*vp);
153 GINAC_ASSERT(is_canonical());
157 // functions overriding virtual functions from bases classes
162 basic * expairseq::duplicate() const
164 debugmsg("expairseq duplicate",LOGLEVEL_DUPLICATE);
165 return new expairseq(*this);
168 void expairseq::print(ostream & os, unsigned upper_precedence) const
170 debugmsg("expairseq print",LOGLEVEL_PRINT);
172 printseq(os,',',precedence,upper_precedence);
176 void expairseq::printraw(ostream & os) const
178 debugmsg("expairseq printraw",LOGLEVEL_PRINT);
181 for (epvector::const_iterator cit=seq.begin(); cit!=seq.end(); ++cit) {
183 (*cit).rest.printraw(os);
185 (*cit).coeff.printraw(os);
191 void expairseq::printtree(ostream & os, unsigned indent) const
193 debugmsg("expairseq printtree",LOGLEVEL_PRINT);
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;
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);
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;
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;
222 for (unsigned i=0; i<hashtabsize; ++i) {
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() << " ";
233 cum_fill += this_bin_fill;
234 cum_fill_sq += this_bin_fill*this_bin_fill;
236 if (this_bin_fill<MAXCOUNT) {
237 ++count[this_bin_fill];
244 double lambda=(1.0*seq.size())/hashtabsize;
245 for (int k=0; k<MAXCOUNT; ++k) {
247 double prob=pow(lambda,k)/fact*exp(-lambda);
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;
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;
257 os << string(indent+delta_indent,' ') << "variance: "
258 << 1.0/hashtabsize*cum_fill_sq-(1.0/hashtabsize*cum_fill)*(1.0/hashtabsize*cum_fill)
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
266 bool expairseq::info(unsigned inf) const
268 return basic::info(inf);
271 int expairseq::nops() const
273 if (overall_coeff.is_equal(default_overall_coeff())) {
279 ex expairseq::op(int const i) const
281 if (unsigned(i)<seq.size()) {
282 return recombine_pair_to_ex(seq[i]);
284 GINAC_ASSERT(!overall_coeff.is_equal(default_overall_coeff()));
285 return overall_coeff;
288 ex & expairseq::let_op(int const i)
290 throw(std::logic_error("let_op not defined for expairseq and derived classes (add,mul,...)"));
293 ex expairseq::eval(int level) const
295 if ((level==1)&&(flags & status_flags::evaluated)) {
299 epvector * vp=evalchildren(level);
304 return (new expairseq(vp,overall_coeff))
305 ->setflag(status_flags::dynallocated |
306 status_flags::evaluated );
309 ex expairseq::evalf(int level) const
311 return thisexpairseq(evalfchildren(level),overall_coeff);
314 ex expairseq::normal(lst &sym_lst, lst &repl_lst, int level) const
316 ex n=thisexpairseq(normalchildren(level),overall_coeff);
317 return n.bp->basic::normal(sym_lst,repl_lst,level);
320 ex expairseq::subs(lst const & ls, lst const & lr) const
322 epvector * vp=subschildren(ls,lr);
326 return thisexpairseq(vp,overall_coeff);
331 int expairseq::compare_same_type(basic const & other) const
333 GINAC_ASSERT(is_of_type(other, expairseq));
334 expairseq const & o=static_cast<expairseq const &>(const_cast<basic &>(other));
338 // compare number of elements
339 if (seq.size() != o.seq.size()) {
340 return (seq.size()<o.seq.size()) ? -1 : 1;
343 // compare overall_coeff
344 cmpval=overall_coeff.compare(o.overall_coeff);
345 if (cmpval!=0) return cmpval;
347 //if (seq.size()==0) return 0; // empty expairseq's are equal
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();
358 for (; (cit1!=last1)&&(cit2!=last2); ++cit1, ++cit2) {
359 cmpval=(*cit1).compare(*cit2);
360 if (cmpval!=0) return cmpval;
363 GINAC_ASSERT(cit1==last1);
364 GINAC_ASSERT(cit2==last2);
367 #ifdef EXPAIRSEQ_USE_HASHTAB
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;
378 // compare individual (sorted) hashtab entries
379 for (unsigned i=0; i<hashtabsize; ++i) {
380 unsigned sz=hashtab[i].size();
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;
396 #endif // def EXPAIRSEQ_USE_HASHTAB
399 bool expairseq::is_equal_same_type(basic const & other) const
401 expairseq const & o=dynamic_cast<expairseq const &>(const_cast<basic &>(other));
403 // compare number of elements
404 if (seq.size() != o.seq.size()) return false;
406 // compare overall_coeff
407 if (!overall_coeff.is_equal(o.overall_coeff)) return false;
409 #ifdef EXPAIRSEQ_USE_HASHTAB
410 // compare number of elements in each hashtab entry
411 if (hashtabsize!=o.hashtabsize) {
412 cout << "this:" << endl;
414 cout << "other:" << endl;
415 other.printtree(cout,0);
418 GINAC_ASSERT(hashtabsize==o.hashtabsize);
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();
426 while (cit1!=last1) {
427 if (!(*cit1).is_equal(*cit2)) return false;
433 #ifdef EXPAIRSEQ_USE_HASHTAB
436 for (unsigned i=0; i<hashtabsize; ++i) {
437 if (hashtab[i].size() != o.hashtab[i].size()) return false;
440 // compare individual sorted hashtab entries
441 for (unsigned i=0; i<hashtabsize; ++i) {
442 unsigned sz=hashtab[i].size();
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;
457 #endif // def EXPAIRSEQ_USE_HASHTAB
460 unsigned expairseq::return_type(void) const
462 return return_types::noncommutative_composite;
465 unsigned expairseq::calchash(void) const
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();
476 v ^= overall_coeff.gethash();
479 // store calculated hash value only if object is already evaluated
480 if (flags & status_flags::evaluated) {
481 setflag(status_flags::hash_calculated);
488 ex expairseq::expand(unsigned options) const
490 epvector * vp=expandchildren(options);
494 return thisexpairseq(vp,overall_coeff);
498 // new virtual functions which can be overridden by derived classes
503 ex expairseq::thisexpairseq(epvector const & v,ex const & oc) const
505 return expairseq(v,oc);
508 ex expairseq::thisexpairseq(epvector * vp, ex const & oc) const
510 return expairseq(vp,oc);
513 void expairseq::printpair(ostream & os, expair const & p, unsigned upper_precedence) const
516 p.rest.bp->print(os,precedence);
518 p.coeff.bp->print(os,precedence);
522 void expairseq::printseq(ostream & os, char delim, unsigned this_precedence,
523 unsigned upper_precedence) const
525 if (this_precedence<=upper_precedence) os << "(";
526 epvector::const_iterator it,it_last;
529 for (it=seq.begin(); it!=it_last; ++it) {
530 printpair(os,*it,this_precedence);
533 printpair(os,*it,this_precedence);
534 if (!overall_coeff.is_equal(default_overall_coeff())) {
535 os << delim << overall_coeff;
537 if (this_precedence<=upper_precedence) os << ")";
540 expair expairseq::split_ex_to_pair(ex const & e) const
542 return expair(e,_ex1());
545 expair expairseq::combine_ex_with_coeff_to_pair(ex const & e,
548 GINAC_ASSERT(is_ex_exactly_of_type(c,numeric));
553 expair expairseq::combine_pair_with_coeff_to_pair(expair const & p,
556 GINAC_ASSERT(is_ex_exactly_of_type(p.coeff,numeric));
557 GINAC_ASSERT(is_ex_exactly_of_type(c,numeric));
559 return expair(p.rest,ex_to_numeric(p.coeff).mul_dyn(ex_to_numeric(c)));
562 ex expairseq::recombine_pair_to_ex(expair const & p) const
564 return lst(p.rest,p.coeff);
567 bool expairseq::expair_needs_further_processing(epp it)
572 ex expairseq::default_overall_coeff(void) const
577 void expairseq::combine_overall_coeff(ex const & c)
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));
584 void expairseq::combine_overall_coeff(ex const & c1, ex const & c2)
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)));
593 bool expairseq::can_make_flat(expair const & p) const
600 // non-virtual functions in this class
603 void expairseq::construct_from_2_ex_via_exvector(ex const & lh, ex const & 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
616 void expairseq::construct_from_2_ex(ex const & lh, ex const & rh)
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);
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
631 #endif // def EXPAIRSEQ_USE_HASHTAB
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);
639 #endif // def EXPAIRSEQ_USE_HASHTAB
640 construct_from_expairseq_ex(ex_to_expairseq(lh),rh);
641 #ifdef EXPAIRSEQ_USE_HASHTAB
643 #endif // def EXPAIRSEQ_USE_HASHTAB
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);
652 #endif // def EXPAIRSEQ_USE_HASHTAB
653 construct_from_expairseq_ex(ex_to_expairseq(rh),lh);
654 #ifdef EXPAIRSEQ_USE_HASHTAB
656 #endif // def EXPAIRSEQ_USE_HASHTAB
660 #ifdef EXPAIRSEQ_USE_HASHTAB
661 if (calc_hashtabsize(2)!=0) {
662 construct_from_2_ex_via_exvector(lh,rh);
666 #endif // def EXPAIRSEQ_USE_HASHTAB
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);
673 combine_overall_coeff(lh);
674 seq.push_back(split_ex_to_pair(rh));
677 if (is_ex_exactly_of_type(rh,numeric)) {
678 combine_overall_coeff(rh);
679 seq.push_back(split_ex_to_pair(lh));
681 expair p1=split_ex_to_pair(lh);
682 expair p2=split_ex_to_pair(rh);
684 int cmpval=p1.rest.compare(p2.rest);
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()
706 void expairseq::construct_from_2_expairseq(expairseq const & s1,
707 expairseq const & s2)
709 combine_overall_coeff(s1.overall_coeff);
710 combine_overall_coeff(s2.overall_coeff);
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();
717 seq.reserve(s1.seq.size()+s2.seq.size());
719 bool needs_further_processing=false;
721 while (first1!=last1 && first2!=last2) {
722 int cmpval=(*first1).rest.compare((*first2).rest);
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;
735 } else if (cmpval<0) {
736 seq.push_back(*first1);
739 seq.push_back(*first2);
744 while (first1!=last1) {
745 seq.push_back(*first1);
748 while (first2!=last2) {
749 seq.push_back(*first2);
753 if (needs_further_processing) {
756 construct_from_epvector(v);
760 void expairseq::construct_from_expairseq_ex(expairseq const & s,
763 combine_overall_coeff(s.overall_coeff);
764 if (is_ex_exactly_of_type(e,numeric)) {
765 combine_overall_coeff(e);
770 epvector::const_iterator first=s.seq.begin();
771 epvector::const_iterator last=s.seq.end();
772 expair p=split_ex_to_pair(e);
774 seq.reserve(s.seq.size()+1);
777 bool needs_further_processing=false;
779 // merge p into s.seq
780 while (first!=last) {
781 int cmpval=(*first).rest.compare(p.rest);
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;
795 } else if (cmpval<0) {
796 seq.push_back(*first);
806 // while loop exited because p was pushed, now push rest of s.seq
807 while (first!=last) {
808 seq.push_back(*first);
812 // while loop exited because s.seq was pushed, now push p
816 if (needs_further_processing) {
819 construct_from_epvector(v);
823 void expairseq::construct_from_exvector(exvector const & v)
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 (+,*) -> (*,^)
831 #ifdef EXPAIRSEQ_USE_HASHTAB
832 combine_same_terms();
835 combine_same_terms_sorted_seq();
836 #endif // def EXPAIRSEQ_USE_HASHTAB
839 void expairseq::construct_from_epvector(epvector const & v)
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 (+,*) -> (*,^)
847 #ifdef EXPAIRSEQ_USE_HASHTAB
848 combine_same_terms();
851 combine_same_terms_sorted_seq();
852 #endif // def EXPAIRSEQ_USE_HASHTAB
857 void expairseq::make_flat(exvector const & v)
859 exvector::const_iterator cit, citend = v.end();
861 // count number of operands which are of same expairseq derived type
862 // and their cumulative number of operands
866 while (cit!=citend) {
867 if (cit->bp->tinfo()==tinfo()) {
869 noperands+=ex_to_expairseq(*cit).seq.size();
874 // reserve seq and coeffseq which will hold all operands
875 seq.reserve(v.size()+noperands-nexpairseqs);
877 // copy elements and split off numerical part
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);
889 if (is_ex_exactly_of_type(*cit,numeric)) {
890 combine_overall_coeff(*cit);
892 seq.push_back(split_ex_to_pair(*cit));
899 cout << "after make flat" << endl;
900 for (epvector::const_iterator cit=seq.begin(); cit!=seq.end(); ++cit) {
901 (*cit).printraw(cout);
908 void expairseq::make_flat(epvector const & v)
910 epvector::const_iterator cit, citend = v.end();
912 // count number of operands which are of same expairseq derived type
913 // and their cumulative number of operands
917 while (cit!=citend) {
918 if (cit->rest.bp->tinfo()==tinfo()) {
920 noperands+=ex_to_expairseq((*cit).rest).seq.size();
925 // reserve seq and coeffseq which will hold all operands
926 seq.reserve(v.size()+noperands-nexpairseqs);
928 // copy elements and split off numerical part
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,
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));
956 epvector * expairseq::bubblesort(epvector::iterator itbegin, epvector::iterator itend)
958 unsigned n=itend-itbegin;
960 epvector * sp=new epvector;
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) {
972 sp->push_back(*last);
976 epvector * expairseq::mergesort(epvector::iterator itbegin, epvector::iterator itend)
978 unsigned n=itend-itbegin;
981 epvector * sp=new epvector;
982 sp->push_back(*itbegin);
986 if (n<16) return bubblesort(itbegin, itend);
989 epvector * s1p=mergesort(itbegin, itbegin+m);
990 epvector * s2p=mergesort(itbegin+m, itend);
992 epvector * sp=new epvector;
993 sp->reserve(s1p->size()+s2p->size());
995 epvector::iterator first1=s1p->begin();
996 epvector::iterator last1=s1p->end();
998 epvector::iterator first2=s2p->begin();
999 epvector::iterator last2=s2p->end();
1001 while (first1 != last1 && first2 != last2) {
1002 if ((*first1).rest.compare((*first2).rest)<0) {
1003 sp->push_back(*first1);
1006 sp->push_back(*first2);
1011 if (first1 != last1) {
1012 while (first1 != last1) {
1013 sp->push_back(*first1);
1017 while (first2 != last2) {
1018 sp->push_back(*first2);
1030 void expairseq::canonicalize(void)
1033 sort(seq.begin(),seq.end(),expair_is_less());
1035 sort(seq.begin(),seq.end(),expair_is_less_old());
1037 if (is_ex_exactly_of_type((*(seq.begin())).rest,numeric)) {
1038 sort(seq.begin(),seq.end(),expair_is_less());
1040 epvector::iterator last_numeric=seq.end();
1043 } while (is_ex_exactly_of_type((*last_numeric).rest,numeric));
1045 sort(last_numeric,seq.end(),expair_is_less());
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) {
1061 cout << "after canonicalize" << endl;
1062 for (epvector::const_iterator cit=seq.begin(); cit!=seq.end(); ++cit) {
1063 (*cit).printraw(cout);
1070 void expairseq::combine_same_terms_sorted_seq(void)
1072 bool needs_further_processing=false;
1074 // combine same terms, drop term with coeff 0
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;
1092 if (!ex_to_numeric((*itin1).coeff).is_zero()) {
1102 if (!ex_to_numeric((*itin1).coeff).is_zero()) {
1109 seq.erase(itout,last);
1114 cout << "after combine" << endl;
1115 for (epvector::const_iterator cit=seq.begin(); cit!=seq.end(); ++cit) {
1116 (*cit).printraw(cout);
1122 if (needs_further_processing) {
1125 construct_from_epvector(v);
1129 #ifdef EXPAIRSEQ_USE_HASHTAB
1131 unsigned expairseq::calc_hashtabsize(unsigned sz) const
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);
1145 unsigned expairseq::calc_hashindex(ex const & e) const
1147 // calculate hashindex
1148 unsigned hash=e.gethash();
1150 if (is_a_numeric_hash(hash)) {
1153 hashindex=hash & hashmask;
1154 // last hashtab entry is reserved for numerics
1155 if (hashindex==hashmask) hashindex=0;
1157 GINAC_ASSERT(hashindex>=0);
1158 GINAC_ASSERT((hashindex<hashtabsize)||(hashtabsize==0));
1162 void expairseq::shrink_hashtab(void)
1164 unsigned new_hashtabsize;
1165 while (hashtabsize!=(new_hashtabsize=calc_hashtabsize(seq.size()))) {
1166 GINAC_ASSERT(new_hashtabsize<hashtabsize);
1167 if (new_hashtabsize==0) {
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());
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;
1188 void expairseq::remove_hashtab_entry(epvector::const_iterator element)
1190 if (hashtabsize==0) return; // nothing to do
1192 // calculate hashindex of element to be deleted
1193 unsigned hashindex=calc_hashindex((*element).rest);
1195 // find it in hashtab and remove it
1196 epplist & eppl=hashtab[hashindex];
1197 epplist::iterator epplit=eppl.begin();
1199 while (epplit!=eppl.end()) {
1200 if (*epplit == element) {
1209 cout << "tried to erase " << element-seq.begin() << endl;
1210 cout << "size " << seq.end()-seq.begin() << endl;
1212 unsigned hashindex=calc_hashindex((*element).rest);
1213 epplist & eppl=hashtab[hashindex];
1214 epplist::iterator epplit=eppl.begin();
1216 while (epplit!=eppl.end()) {
1217 if (*epplit == element) {
1224 GINAC_ASSERT(erased);
1226 GINAC_ASSERT(erased);
1229 void expairseq::move_hashtab_entry(epvector::const_iterator oldpos,
1230 epvector::iterator newpos)
1232 GINAC_ASSERT(hashtabsize!=0);
1234 // calculate hashindex of element which was moved
1235 unsigned hashindex=calc_hashindex((*newpos).rest);
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) {
1247 GINAC_ASSERT(epplit!=eppl.end());
1250 void expairseq::sorted_insert(epplist & eppl, epp elem)
1252 epplist::iterator current=eppl.begin();
1253 while ((current!=eppl.end())&&((*(*current)).is_less(*elem))) {
1256 eppl.insert(current,elem);
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)
1264 epp current=seq.begin();
1266 while (current!=first_numeric) {
1267 if (is_ex_exactly_of_type((*current).rest,numeric)) {
1269 iter_swap(current,first_numeric);
1271 // calculate hashindex
1272 unsigned currenthashindex=calc_hashindex((*current).rest);
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;
1281 if (epplit==eppl.end()) {
1282 // no matching expair found, append this to end of list
1283 sorted_insert(eppl,current);
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));
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);
1294 if (first_numeric!=last_non_zero) iter_swap(first_numeric,current);
1297 // test if combined term has coeff 0 and can be removed is done later
1298 touched[(*epplit)-seq.begin()]=true;
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)
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();
1313 while (current!=first_numeric) {
1317 } else if (!ex_to_numeric((*current).coeff).is_equal(_num0())) {
1321 remove_hashtab_entry(current);
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);
1327 bool numeric_swapped=first_numeric!=last_non_zero;
1328 if (numeric_swapped) iter_swap(first_numeric,current);
1329 epvector::iterator changed_entry;
1331 if (numeric_swapped) {
1332 changed_entry=first_numeric;
1334 changed_entry=last_non_zero;
1340 if (first_numeric!=current) {
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()];
1353 GINAC_ASSERT(i==current-seq.begin());
1356 bool expairseq::has_coeff_0(void) const
1358 for (epvector::const_iterator cit=seq.begin(); cit!=seq.end(); ++cit) {
1359 if ((*cit).coeff.is_equal(_ex0())) {
1366 void expairseq::add_numerics_to_hashtab(epvector::iterator first_numeric,
1367 epvector::const_iterator last_non_zero)
1369 if (first_numeric==seq.end()) return; // no numerics
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);
1379 void expairseq::combine_same_terms(void)
1381 // combine same terms, drop term with coeff 0, move numerics to end
1383 // calculate size of hashtab
1384 hashtabsize=calc_hashtabsize(seq.size());
1386 // hashtabsize is a power of 2
1387 hashmask=hashtabsize-1;
1391 hashtab.resize(hashtabsize);
1393 if (hashtabsize==0) {
1395 combine_same_terms_sorted_seq();
1396 GINAC_ASSERT(!has_coeff_0());
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;
1405 vector<bool> touched;
1406 touched.reserve(seq.size());
1407 for (unsigned i=0; i<seq.size(); ++i) touched[i]=false;
1409 unsigned number_of_zeroes=0;
1411 GINAC_ASSERT(!has_coeff_0());
1412 build_hashtab_and_combine(first_numeric,last_non_zero,touched,number_of_zeroes);
1414 cout << "in combine:" << endl;
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;
1422 cout << "end in combine" << endl;
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);
1430 cout << "in combine after drop:" << endl;
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;
1438 cout << "end in combine after drop" << endl;
1442 add_numerics_to_hashtab(first_numeric,last_non_zero);
1444 // pop zero elements
1445 for (unsigned i=0; i<number_of_zeroes; ++i) {
1449 // shrink hashtabsize to calculated value
1450 GINAC_ASSERT(!has_coeff_0());
1454 GINAC_ASSERT(!has_coeff_0());
1457 #endif // def EXPAIRSEQ_USE_HASHTAB
1459 bool expairseq::is_canonical() const
1461 if (seq.size()<=1) return 1;
1463 #ifdef EXPAIRSEQ_USE_HASHTAB
1464 if (hashtabsize>0) return 1; // not canoncalized
1465 #endif // def EXPAIRSEQ_USE_HASHTAB
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);
1478 printpair(cout,*it,0);
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);
1494 epvector * expairseq::expandchildren(unsigned options) const
1496 epvector::const_iterator last=seq.end();
1497 epvector::const_iterator cit=seq.begin();
1499 ex const & expanded_ex=(*cit).rest.expand(options);
1500 if (!are_ex_trivially_equal((*cit).rest,expanded_ex)) {
1502 // something changed, copy seq, eval and return it
1503 epvector *s=new epvector;
1504 s->reserve(seq.size());
1506 // copy parts of seq which are known not to have changed
1507 epvector::const_iterator cit2=seq.begin();
1509 s->push_back(*cit2);
1512 // copy first changed element
1513 s->push_back(combine_ex_with_coeff_to_pair(expanded_ex,
1517 while (cit2!=last) {
1518 s->push_back(combine_ex_with_coeff_to_pair((*cit2).rest.expand(options),
1527 return 0; // nothing has changed
1530 epvector * expairseq::evalchildren(int level) const
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)
1539 if (level == -max_recursion_level) {
1540 throw(std::runtime_error("max recursion level reached"));
1544 epvector::const_iterator last=seq.end();
1545 epvector::const_iterator cit=seq.begin();
1547 ex const & evaled_ex=(*cit).rest.eval(level);
1548 if (!are_ex_trivially_equal((*cit).rest,evaled_ex)) {
1550 // something changed, copy seq, eval and return it
1551 epvector *s=new epvector;
1552 s->reserve(seq.size());
1554 // copy parts of seq which are known not to have changed
1555 epvector::const_iterator cit2=seq.begin();
1557 s->push_back(*cit2);
1560 // copy first changed element
1561 s->push_back(combine_ex_with_coeff_to_pair(evaled_ex,
1565 while (cit2!=last) {
1566 s->push_back(combine_ex_with_coeff_to_pair((*cit2).rest.eval(level),
1575 return 0; // nothing has changed
1578 epvector expairseq::evalfchildren(int level) const
1581 s.reserve(seq.size());
1586 if (level == -max_recursion_level) {
1587 throw(std::runtime_error("max recursion level reached"));
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),
1597 epvector expairseq::normalchildren(int level) const
1600 s.reserve(seq.size());
1605 if (level == -max_recursion_level) {
1606 throw(std::runtime_error("max recursion level reached"));
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),
1616 epvector expairseq::diffchildren(symbol const & y) const
1619 s.reserve(seq.size());
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),
1628 epvector * expairseq::subschildren(lst const & ls, lst const & lr) const
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());
1635 epvector::const_iterator last=seq.end();
1636 epvector::const_iterator cit=seq.begin();
1638 ex const & subsed_ex=(*cit).rest.subs(ls,lr);
1639 if (!are_ex_trivially_equal((*cit).rest,subsed_ex)) {
1641 // something changed, copy seq, subs and return it
1642 epvector *s=new epvector;
1643 s->reserve(seq.size());
1645 // copy parts of seq which are known not to have changed
1646 epvector::const_iterator cit2=seq.begin();
1648 s->push_back(*cit2);
1651 // copy first changed element
1652 s->push_back(combine_ex_with_coeff_to_pair(subsed_ex,
1656 while (cit2!=last) {
1657 s->push_back(combine_ex_with_coeff_to_pair((*cit2).rest.subs(ls,lr),
1666 return 0; // nothing has changed
1670 // static member variables
1675 unsigned expairseq::precedence=10;
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
1687 const expairseq some_expairseq;
1688 type_info const & typeid_expairseq=typeid(some_expairseq);
1690 #ifndef NO_GINAC_NAMESPACE
1691 } // namespace GiNaC
1692 #endif // ndef NO_GINAC_NAMESPACE