1 /** @file expairseq.cpp
3 * Implementation of sequences of expression pairs. */
6 * GiNaC Copyright (C) 1999-2001 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
27 #include "expairseq.h"
33 #if EXPAIRSEQ_USE_HASHTAB
35 #endif // EXPAIRSEQ_USE_HASHTAB
39 GINAC_IMPLEMENT_REGISTERED_CLASS_NO_CTORS(expairseq, basic)
48 bool operator()(const epp &lh, const epp &rh) const
50 return (*lh).is_less(*rh);
55 // default ctor, dtor, copy ctor assignment operator and helpers
60 expairseq::expairseq(const expairseq &other)
62 debugmsg("expairseq copy ctor",LOGLEVEL_CONSTRUCT);
66 const expairseq &expairseq::operator=(const expairseq &other)
68 debugmsg("expairseq operator=",LOGLEVEL_ASSIGNMENT);
78 /** For use by copy ctor and assignment operator. */
79 void expairseq::copy(const expairseq &other)
81 inherited::copy(other);
83 overall_coeff = other.overall_coeff;
84 #if EXPAIRSEQ_USE_HASHTAB
86 hashtabsize = other.hashtabsize;
88 hashmask = other.hashmask;
89 hashtab.resize(hashtabsize);
90 epvector::const_iterator osb = other.seq.begin();
91 for (unsigned i=0; i<hashtabsize; ++i) {
93 for (epplist::const_iterator cit=other.hashtab[i].begin();
94 cit!=other.hashtab[i].end(); ++cit) {
95 hashtab[i].push_back(seq.begin()+((*cit)-osb));
101 #endif // EXPAIRSEQ_USE_HASHTAB
104 void expairseq::destroy(bool call_parent)
107 basic::destroy(call_parent);
114 expairseq::expairseq(const ex &lh, const ex &rh) : inherited(TINFO_expairseq)
116 debugmsg("expairseq ctor from ex,ex",LOGLEVEL_CONSTRUCT);
117 construct_from_2_ex(lh,rh);
118 GINAC_ASSERT(is_canonical());
121 expairseq::expairseq(const exvector &v) : inherited(TINFO_expairseq)
123 debugmsg("expairseq ctor from exvector",LOGLEVEL_CONSTRUCT);
124 construct_from_exvector(v);
125 GINAC_ASSERT(is_canonical());
128 expairseq::expairseq(const epvector &v, const ex &oc)
129 : inherited(TINFO_expairseq), overall_coeff(oc)
131 debugmsg("expairseq ctor from epvector,ex",LOGLEVEL_CONSTRUCT);
132 construct_from_epvector(v);
133 GINAC_ASSERT(is_canonical());
136 expairseq::expairseq(epvector *vp, const ex &oc)
137 : inherited(TINFO_expairseq), overall_coeff(oc)
139 debugmsg("expairseq ctor from epvector *,ex",LOGLEVEL_CONSTRUCT);
141 construct_from_epvector(*vp);
143 GINAC_ASSERT(is_canonical());
150 /** Construct object from archive_node. */
151 expairseq::expairseq(const archive_node &n, const lst &sym_lst) : inherited(n, sym_lst)
152 #if EXPAIRSEQ_USE_HASHTAB
156 debugmsg("expairseq ctor from archive_node", LOGLEVEL_CONSTRUCT);
157 for (unsigned int i=0; true; i++) {
160 if (n.find_ex("rest", rest, sym_lst, i) && n.find_ex("coeff", coeff, sym_lst, i))
161 seq.push_back(expair(rest, coeff));
165 n.find_ex("overall_coeff", overall_coeff, sym_lst);
168 /** Unarchive the object. */
169 ex expairseq::unarchive(const archive_node &n, const lst &sym_lst)
171 return (new expairseq(n, sym_lst))->setflag(status_flags::dynallocated);
174 /** Archive the object. */
175 void expairseq::archive(archive_node &n) const
177 inherited::archive(n);
178 epvector::const_iterator i = seq.begin(), iend = seq.end();
180 n.add_ex("rest", i->rest);
181 n.add_ex("coeff", i->coeff);
184 n.add_ex("overall_coeff", overall_coeff);
188 // functions overriding virtual functions from bases classes
193 basic *expairseq::duplicate() const
195 debugmsg("expairseq duplicate",LOGLEVEL_DUPLICATE);
196 return new expairseq(*this);
199 void expairseq::print(std::ostream &os, unsigned upper_precedence) const
201 debugmsg("expairseq print",LOGLEVEL_PRINT);
203 printseq(os,',',precedence,upper_precedence);
207 void expairseq::printraw(std::ostream &os) const
209 debugmsg("expairseq printraw",LOGLEVEL_PRINT);
210 os << class_name() << "(";
211 for (epvector::const_iterator cit=seq.begin(); cit!=seq.end(); ++cit) {
213 (*cit).rest.printraw(os);
215 (*cit).coeff.printraw(os);
221 void expairseq::printtree(std::ostream &os, unsigned indent) const
223 debugmsg("expairseq printtree",LOGLEVEL_PRINT);
225 os << std::string(indent,' ') << "type=" << class_name()
226 << ", hash=" << hashvalue
227 << " (0x" << std::hex << hashvalue << std::dec << ")"
228 << ", flags=" << flags
229 << ", nops=" << nops() << std::endl;
230 for (unsigned i=0; i<seq.size(); ++i) {
231 seq[i].rest.printtree(os,indent+delta_indent);
232 seq[i].coeff.printtree(os,indent+delta_indent);
234 os << std::string(indent+delta_indent,' ') << "-----" << std::endl;
236 if (!overall_coeff.is_equal(default_overall_coeff())) {
237 os << std::string(indent+delta_indent,' ') << "-----" << std::endl;
238 os << std::string(indent+delta_indent,' ') << "overall_coeff" << std::endl;
239 overall_coeff.printtree(os,indent+delta_indent);
241 os << std::string(indent+delta_indent,' ') << "=====" << std::endl;
242 #if EXPAIRSEQ_USE_HASHTAB
243 os << std::string(indent+delta_indent,' ')
244 << "hashtab size " << hashtabsize << std::endl;
245 if (hashtabsize==0) return;
247 unsigned count[MAXCOUNT+1];
248 for (int i=0; i<MAXCOUNT+1; ++i)
250 unsigned this_bin_fill;
251 unsigned cum_fill_sq = 0;
252 unsigned cum_fill = 0;
253 for (unsigned i=0; i<hashtabsize; ++i) {
255 if (hashtab[i].size()>0) {
256 os << std::string(indent+delta_indent,' ')
257 << "bin " << i << " with entries ";
258 for (epplist::const_iterator it=hashtab[i].begin();
259 it!=hashtab[i].end(); ++it) {
260 os << *it-seq.begin() << " ";
264 cum_fill += this_bin_fill;
265 cum_fill_sq += this_bin_fill*this_bin_fill;
267 if (this_bin_fill<MAXCOUNT)
268 ++count[this_bin_fill];
274 double lambda = (1.0*seq.size())/hashtabsize;
275 for (int k=0; k<MAXCOUNT; ++k) {
278 double prob = std::pow(lambda,k)/fact * std::exp(-lambda);
280 os << std::string(indent+delta_indent,' ') << "bins with " << k << " entries: "
281 << int(1000.0*count[k]/hashtabsize)/10.0 << "% (expected: "
282 << int(prob*1000)/10.0 << ")" << std::endl;
284 os << std::string(indent+delta_indent,' ') << "bins with more entries: "
285 << int(1000.0*count[MAXCOUNT]/hashtabsize)/10.0 << "% (expected: "
286 << int((1-cum_prob)*1000)/10.0 << ")" << std::endl;
288 os << std::string(indent+delta_indent,' ') << "variance: "
289 << 1.0/hashtabsize*cum_fill_sq-(1.0/hashtabsize*cum_fill)*(1.0/hashtabsize*cum_fill)
291 os << std::string(indent+delta_indent,' ') << "average fill: "
292 << (1.0*cum_fill)/hashtabsize
293 << " (should be equal to " << (1.0*seq.size())/hashtabsize << ")" << std::endl;
294 #endif // EXPAIRSEQ_USE_HASHTAB
297 bool expairseq::info(unsigned inf) const
299 return inherited::info(inf);
302 unsigned expairseq::nops() const
304 if (overall_coeff.is_equal(default_overall_coeff()))
310 ex expairseq::op(int i) const
312 if (unsigned(i)<seq.size())
313 return recombine_pair_to_ex(seq[i]);
314 GINAC_ASSERT(!overall_coeff.is_equal(default_overall_coeff()));
315 return overall_coeff;
318 ex &expairseq::let_op(int i)
320 throw(std::logic_error("let_op not defined for expairseq and derived classes (add,mul,...)"));
323 ex expairseq::eval(int level) const
325 if ((level==1) && (flags &status_flags::evaluated))
328 epvector *vp = evalchildren(level);
332 return (new expairseq(vp,overall_coeff))->setflag(status_flags::dynallocated | status_flags::evaluated);
335 ex expairseq::evalf(int level) const
337 return thisexpairseq(evalfchildren(level),overall_coeff.evalf(level-1));
340 ex expairseq::normal(lst &sym_lst, lst &repl_lst, int level) const
342 ex n = thisexpairseq(normalchildren(level),overall_coeff);
343 return n.bp->basic::normal(sym_lst,repl_lst,level);
346 ex expairseq::subs(const lst &ls, const lst &lr) const
348 epvector *vp = subschildren(ls,lr);
352 return thisexpairseq(vp,overall_coeff);
357 /** Implementation of ex::diff() for an expairseq.
358 * It differentiates all elements of the sequence.
360 ex expairseq::derivative(const symbol &s) const
362 return thisexpairseq(diffchildren(s),overall_coeff);
365 int expairseq::compare_same_type(const basic &other) const
367 GINAC_ASSERT(is_of_type(other, expairseq));
368 const expairseq &o = static_cast<const expairseq &>(const_cast<basic &>(other));
372 // compare number of elements
373 if (seq.size() != o.seq.size())
374 return (seq.size()<o.seq.size()) ? -1 : 1;
376 // compare overall_coeff
377 cmpval = overall_coeff.compare(o.overall_coeff);
381 #if EXPAIRSEQ_USE_HASHTAB
382 GINAC_ASSERT(hashtabsize==o.hashtabsize);
383 if (hashtabsize==0) {
384 #endif // EXPAIRSEQ_USE_HASHTAB
385 epvector::const_iterator cit1 = seq.begin();
386 epvector::const_iterator cit2 = o.seq.begin();
387 epvector::const_iterator last1 = seq.end();
388 epvector::const_iterator last2 = o.seq.end();
390 for (; (cit1!=last1)&&(cit2!=last2); ++cit1, ++cit2) {
391 cmpval = (*cit1).compare(*cit2);
392 if (cmpval!=0) return cmpval;
395 GINAC_ASSERT(cit1==last1);
396 GINAC_ASSERT(cit2==last2);
399 #if EXPAIRSEQ_USE_HASHTAB
402 // compare number of elements in each hashtab entry
403 for (unsigned i=0; i<hashtabsize; ++i) {
404 unsigned cursize=hashtab[i].size();
405 if (cursize != o.hashtab[i].size())
406 return (cursize < o.hashtab[i].size()) ? -1 : 1;
409 // compare individual (sorted) hashtab entries
410 for (unsigned i=0; i<hashtabsize; ++i) {
411 unsigned sz = hashtab[i].size();
413 const epplist &eppl1 = hashtab[i];
414 const epplist &eppl2 = o.hashtab[i];
415 epplist::const_iterator it1 = eppl1.begin();
416 epplist::const_iterator it2 = eppl2.begin();
417 while (it1!=eppl1.end()) {
418 cmpval = (*(*it1)).compare(*(*it2));
428 #endif // EXPAIRSEQ_USE_HASHTAB
431 bool expairseq::is_equal_same_type(const basic &other) const
433 const expairseq &o = dynamic_cast<const expairseq &>(const_cast<basic &>(other));
435 // compare number of elements
436 if (seq.size()!=o.seq.size())
439 // compare overall_coeff
440 if (!overall_coeff.is_equal(o.overall_coeff))
443 #if EXPAIRSEQ_USE_HASHTAB
444 // compare number of elements in each hashtab entry
445 if (hashtabsize!=o.hashtabsize) {
446 std::cout << "this:" << std::endl;
447 printtree(std::cout,0);
448 std::cout << "other:" << std::endl;
449 other.printtree(std::cout,0);
452 GINAC_ASSERT(hashtabsize==o.hashtabsize);
454 if (hashtabsize==0) {
455 #endif // EXPAIRSEQ_USE_HASHTAB
456 epvector::const_iterator cit1 = seq.begin();
457 epvector::const_iterator cit2 = o.seq.begin();
458 epvector::const_iterator last1 = seq.end();
460 while (cit1!=last1) {
461 if (!(*cit1).is_equal(*cit2)) return false;
467 #if EXPAIRSEQ_USE_HASHTAB
470 for (unsigned i=0; i<hashtabsize; ++i) {
471 if (hashtab[i].size() != o.hashtab[i].size())
475 // compare individual sorted hashtab entries
476 for (unsigned i=0; i<hashtabsize; ++i) {
477 unsigned sz = hashtab[i].size();
479 const epplist &eppl1 = hashtab[i];
480 const epplist &eppl2 = o.hashtab[i];
481 epplist::const_iterator it1 = eppl1.begin();
482 epplist::const_iterator it2 = eppl2.begin();
483 while (it1!=eppl1.end()) {
484 if (!(*(*it1)).is_equal(*(*it2))) return false;
492 #endif // EXPAIRSEQ_USE_HASHTAB
495 unsigned expairseq::return_type(void) const
497 return return_types::noncommutative_composite;
500 unsigned expairseq::calchash(void) const
502 unsigned v = golden_ratio_hash(tinfo());
503 for (epvector::const_iterator cit=seq.begin(); cit!=seq.end(); ++cit) {
504 #if !EXPAIRSEQ_USE_HASHTAB
505 v = rotate_left_31(v); // rotation would spoil commutativity
506 #endif // EXPAIRSEQ_USE_HASHTAB
507 v ^= cit->rest.gethash();
508 #if !EXPAIRSEQ_USE_HASHTAB
509 v = rotate_left_31(v);
510 v ^= cit->coeff.gethash();
511 #endif // EXPAIRSEQ_USE_HASHTAB
514 v ^= overall_coeff.gethash();
517 // store calculated hash value only if object is already evaluated
518 if (flags &status_flags::evaluated) {
519 setflag(status_flags::hash_calculated);
526 ex expairseq::expand(unsigned options) const
528 epvector *vp = expandchildren(options);
530 // the terms have not changed, so it is safe to declare this expanded
531 setflag(status_flags::expanded);
535 return thisexpairseq(vp,overall_coeff);
539 // new virtual functions which can be overridden by derived classes
544 /** Create an object of this type.
545 * This method works similar to a constructor. It is useful because expairseq
546 * has (at least) two possible different semantics but we want to inherit
547 * methods thus avoiding code duplication. Sometimes a method in expairseq
548 * has to create a new one of the same semantics, which cannot be done by a
549 * ctor because the name (add, mul,...) is unknown on the expaiseq level. In
550 * order for this trick to work a derived class must of course override this
552 ex expairseq::thisexpairseq(const epvector &v, const ex &oc) const
554 return expairseq(v,oc);
557 ex expairseq::thisexpairseq(epvector *vp, const ex &oc) const
559 return expairseq(vp,oc);
562 void expairseq::printpair(std::ostream &os, const expair &p, unsigned upper_precedence) const
565 p.rest.bp->print(os,precedence);
567 p.coeff.bp->print(os,precedence);
571 void expairseq::printseq(std::ostream &os, char delim,
572 unsigned this_precedence,
573 unsigned upper_precedence) const
575 if (this_precedence<=upper_precedence)
577 epvector::const_iterator it,it_last;
580 for (it=seq.begin(); it!=it_last; ++it) {
581 printpair(os,*it,this_precedence);
584 printpair(os,*it,this_precedence);
585 if (!overall_coeff.is_equal(default_overall_coeff()))
586 os << delim << overall_coeff;
588 if (this_precedence<=upper_precedence)
593 /** Form an expair from an ex, using the corresponding semantics.
594 * @see expairseq::recombine_pair_to_ex() */
595 expair expairseq::split_ex_to_pair(const ex &e) const
597 return expair(e,_ex1());
601 expair expairseq::combine_ex_with_coeff_to_pair(const ex &e,
604 GINAC_ASSERT(is_ex_exactly_of_type(c,numeric));
610 expair expairseq::combine_pair_with_coeff_to_pair(const expair &p,
613 GINAC_ASSERT(is_ex_exactly_of_type(p.coeff,numeric));
614 GINAC_ASSERT(is_ex_exactly_of_type(c,numeric));
616 return expair(p.rest,ex_to_numeric(p.coeff).mul_dyn(ex_to_numeric(c)));
620 /** Form an ex out of an expair, using the corresponding semantics.
621 * @see expairseq::split_ex_to_pair() */
622 ex expairseq::recombine_pair_to_ex(const expair &p) const
624 return lst(p.rest,p.coeff);
627 bool expairseq::expair_needs_further_processing(epp it)
629 #if EXPAIRSEQ_USE_HASHTAB
630 //# error "FIXME: expair_needs_further_processing not yet implemented for hashtabs, sorry. A.F."
631 #endif // EXPAIRSEQ_USE_HASHTAB
635 ex expairseq::default_overall_coeff(void) const
640 void expairseq::combine_overall_coeff(const ex &c)
642 GINAC_ASSERT(is_ex_exactly_of_type(overall_coeff,numeric));
643 GINAC_ASSERT(is_ex_exactly_of_type(c,numeric));
644 overall_coeff = ex_to_numeric(overall_coeff).add_dyn(ex_to_numeric(c));
647 void expairseq::combine_overall_coeff(const ex &c1, const ex &c2)
649 GINAC_ASSERT(is_ex_exactly_of_type(overall_coeff,numeric));
650 GINAC_ASSERT(is_ex_exactly_of_type(c1,numeric));
651 GINAC_ASSERT(is_ex_exactly_of_type(c2,numeric));
652 overall_coeff = ex_to_numeric(overall_coeff).
653 add_dyn(ex_to_numeric(c1).mul(ex_to_numeric(c2)));
656 bool expairseq::can_make_flat(const expair &p) const
663 // non-virtual functions in this class
666 void expairseq::construct_from_2_ex_via_exvector(const ex &lh, const ex &rh)
672 construct_from_exvector(v);
673 #if EXPAIRSEQ_USE_HASHTAB
674 GINAC_ASSERT((hashtabsize==0)||(hashtabsize>=minhashtabsize));
675 GINAC_ASSERT(hashtabsize==calc_hashtabsize(seq.size()));
676 #endif // EXPAIRSEQ_USE_HASHTAB
679 void expairseq::construct_from_2_ex(const ex &lh, const ex &rh)
681 if (lh.bp->tinfo()==tinfo()) {
682 if (rh.bp->tinfo()==tinfo()) {
683 #if EXPAIRSEQ_USE_HASHTAB
684 unsigned totalsize = ex_to_expairseq(lh).seq.size() +
685 ex_to_expairseq(rh).seq.size();
686 if (calc_hashtabsize(totalsize)!=0) {
687 construct_from_2_ex_via_exvector(lh,rh);
689 #endif // EXPAIRSEQ_USE_HASHTAB
690 construct_from_2_expairseq(ex_to_expairseq(lh),
691 ex_to_expairseq(rh));
692 #if EXPAIRSEQ_USE_HASHTAB
694 #endif // EXPAIRSEQ_USE_HASHTAB
697 #if EXPAIRSEQ_USE_HASHTAB
698 unsigned totalsize = ex_to_expairseq(lh).seq.size()+1;
699 if (calc_hashtabsize(totalsize)!=0) {
700 construct_from_2_ex_via_exvector(lh, rh);
702 #endif // EXPAIRSEQ_USE_HASHTAB
703 construct_from_expairseq_ex(ex_to_expairseq(lh), rh);
704 #if EXPAIRSEQ_USE_HASHTAB
706 #endif // EXPAIRSEQ_USE_HASHTAB
709 } else if (rh.bp->tinfo()==tinfo()) {
710 #if EXPAIRSEQ_USE_HASHTAB
711 unsigned totalsize=ex_to_expairseq(rh).seq.size()+1;
712 if (calc_hashtabsize(totalsize)!=0) {
713 construct_from_2_ex_via_exvector(lh,rh);
715 #endif // EXPAIRSEQ_USE_HASHTAB
716 construct_from_expairseq_ex(ex_to_expairseq(rh),lh);
717 #if EXPAIRSEQ_USE_HASHTAB
719 #endif // EXPAIRSEQ_USE_HASHTAB
723 #if EXPAIRSEQ_USE_HASHTAB
724 if (calc_hashtabsize(2)!=0) {
725 construct_from_2_ex_via_exvector(lh,rh);
729 #endif // EXPAIRSEQ_USE_HASHTAB
731 if (is_ex_exactly_of_type(lh,numeric)) {
732 if (is_ex_exactly_of_type(rh,numeric)) {
733 combine_overall_coeff(lh);
734 combine_overall_coeff(rh);
736 combine_overall_coeff(lh);
737 seq.push_back(split_ex_to_pair(rh));
740 if (is_ex_exactly_of_type(rh,numeric)) {
741 combine_overall_coeff(rh);
742 seq.push_back(split_ex_to_pair(lh));
744 expair p1 = split_ex_to_pair(lh);
745 expair p2 = split_ex_to_pair(rh);
747 int cmpval = p1.rest.compare(p2.rest);
749 p1.coeff=ex_to_numeric(p1.coeff).add_dyn(ex_to_numeric(p2.coeff));
750 if (!ex_to_numeric(p1.coeff).is_zero()) {
751 // no further processing is necessary, since this
752 // one element will usually be recombined in eval()
769 void expairseq::construct_from_2_expairseq(const expairseq &s1,
772 combine_overall_coeff(s1.overall_coeff);
773 combine_overall_coeff(s2.overall_coeff);
775 epvector::const_iterator first1 = s1.seq.begin();
776 epvector::const_iterator last1 = s1.seq.end();
777 epvector::const_iterator first2 = s2.seq.begin();
778 epvector::const_iterator last2 = s2.seq.end();
780 seq.reserve(s1.seq.size()+s2.seq.size());
782 bool needs_further_processing=false;
784 while (first1!=last1 && first2!=last2) {
785 int cmpval = (*first1).rest.compare((*first2).rest);
788 const numeric &newcoeff = ex_to_numeric((*first1).coeff).
789 add(ex_to_numeric((*first2).coeff));
790 if (!newcoeff.is_zero()) {
791 seq.push_back(expair((*first1).rest,newcoeff));
792 if (expair_needs_further_processing(seq.end()-1)) {
793 needs_further_processing = true;
798 } else if (cmpval<0) {
799 seq.push_back(*first1);
802 seq.push_back(*first2);
807 while (first1!=last1) {
808 seq.push_back(*first1);
811 while (first2!=last2) {
812 seq.push_back(*first2);
816 if (needs_further_processing) {
819 construct_from_epvector(v);
823 void expairseq::construct_from_expairseq_ex(const expairseq &s,
826 combine_overall_coeff(s.overall_coeff);
827 if (is_ex_exactly_of_type(e,numeric)) {
828 combine_overall_coeff(e);
833 epvector::const_iterator first = s.seq.begin();
834 epvector::const_iterator last = s.seq.end();
835 expair p = split_ex_to_pair(e);
837 seq.reserve(s.seq.size()+1);
838 bool p_pushed = false;
840 bool needs_further_processing=false;
842 // merge p into s.seq
843 while (first!=last) {
844 int cmpval=(*first).rest.compare(p.rest);
847 const numeric &newcoeff = ex_to_numeric((*first).coeff).
848 add(ex_to_numeric(p.coeff));
849 if (!newcoeff.is_zero()) {
850 seq.push_back(expair((*first).rest,newcoeff));
851 if (expair_needs_further_processing(seq.end()-1)) {
852 needs_further_processing = true;
858 } else if (cmpval<0) {
859 seq.push_back(*first);
869 // while loop exited because p was pushed, now push rest of s.seq
870 while (first!=last) {
871 seq.push_back(*first);
875 // while loop exited because s.seq was pushed, now push p
879 if (needs_further_processing) {
882 construct_from_epvector(v);
886 void expairseq::construct_from_exvector(const exvector &v)
888 // simplifications: +(a,+(b,c),d) -> +(a,b,c,d) (associativity)
889 // +(d,b,c,a) -> +(a,b,c,d) (canonicalization)
890 // +(...,x,*(x,c1),*(x,c2)) -> +(...,*(x,1+c1+c2)) (c1, c2 numeric())
891 // (same for (+,*) -> (*,^)
894 #if EXPAIRSEQ_USE_HASHTAB
895 combine_same_terms();
898 combine_same_terms_sorted_seq();
899 #endif // EXPAIRSEQ_USE_HASHTAB
903 void expairseq::construct_from_epvector(const epvector &v)
905 // simplifications: +(a,+(b,c),d) -> +(a,b,c,d) (associativity)
906 // +(d,b,c,a) -> +(a,b,c,d) (canonicalization)
907 // +(...,x,*(x,c1),*(x,c2)) -> +(...,*(x,1+c1+c2)) (c1, c2 numeric())
908 // (same for (+,*) -> (*,^)
911 #if EXPAIRSEQ_USE_HASHTAB
912 combine_same_terms();
915 combine_same_terms_sorted_seq();
916 #endif // EXPAIRSEQ_USE_HASHTAB
920 /** Combine this expairseq with argument exvector.
921 * It cares for associativity as well as for special handling of numerics. */
922 void expairseq::make_flat(const exvector &v)
924 exvector::const_iterator cit;
926 // count number of operands which are of same expairseq derived type
927 // and their cumulative number of operands
932 while (cit!=v.end()) {
933 if (cit->bp->tinfo()==this->tinfo()) {
935 noperands += ex_to_expairseq(*cit).seq.size();
940 // reserve seq and coeffseq which will hold all operands
941 seq.reserve(v.size()+noperands-nexpairseqs);
943 // copy elements and split off numerical part
945 while (cit!=v.end()) {
946 if (cit->bp->tinfo()==this->tinfo()) {
947 const expairseq &subseqref = ex_to_expairseq(*cit);
948 combine_overall_coeff(subseqref.overall_coeff);
949 epvector::const_iterator cit_s = subseqref.seq.begin();
950 while (cit_s!=subseqref.seq.end()) {
951 seq.push_back(*cit_s);
955 if (is_ex_exactly_of_type(*cit,numeric))
956 combine_overall_coeff(*cit);
958 seq.push_back(split_ex_to_pair(*cit));
966 /** Combine this expairseq with argument epvector.
967 * It cares for associativity as well as for special handling of numerics. */
968 void expairseq::make_flat(const epvector &v)
970 epvector::const_iterator cit;
972 // count number of operands which are of same expairseq derived type
973 // and their cumulative number of operands
978 while (cit!=v.end()) {
979 if (cit->rest.bp->tinfo()==this->tinfo()) {
981 noperands += ex_to_expairseq((*cit).rest).seq.size();
986 // reserve seq and coeffseq which will hold all operands
987 seq.reserve(v.size()+noperands-nexpairseqs);
989 // copy elements and split off numerical part
991 while (cit!=v.end()) {
992 if (cit->rest.bp->tinfo()==this->tinfo() &&
993 this->can_make_flat(*cit)) {
994 const expairseq &subseqref = ex_to_expairseq((*cit).rest);
995 combine_overall_coeff(ex_to_numeric(subseqref.overall_coeff),
996 ex_to_numeric((*cit).coeff));
997 epvector::const_iterator cit_s = subseqref.seq.begin();
998 while (cit_s!=subseqref.seq.end()) {
999 seq.push_back(expair((*cit_s).rest,
1000 ex_to_numeric((*cit_s).coeff).mul_dyn(ex_to_numeric((*cit).coeff))));
1001 //seq.push_back(combine_pair_with_coeff_to_pair(*cit_s,
1006 if (cit->is_canonical_numeric())
1007 combine_overall_coeff(cit->rest);
1009 seq.push_back(*cit);
1017 /** Brings this expairseq into a sorted (canonical) form. */
1018 void expairseq::canonicalize(void)
1021 sort(seq.begin(),seq.end(),expair_is_less());
1026 /** Compact a presorted expairseq by combining all matching expairs to one
1027 * each. On an add object, this is responsible for 2*x+3*x+y -> 5*x+y, for
1029 void expairseq::combine_same_terms_sorted_seq(void)
1031 bool needs_further_processing = false;
1034 epvector::iterator itin1 = seq.begin();
1035 epvector::iterator itin2 = itin1+1;
1036 epvector::iterator itout = itin1;
1037 epvector::iterator last = seq.end();
1038 // must_copy will be set to true the first time some combination is
1039 // possible from then on the sequence has changed and must be compacted
1040 bool must_copy = false;
1041 while (itin2!=last) {
1042 if ((*itin1).rest.compare((*itin2).rest)==0) {
1043 (*itin1).coeff = ex_to_numeric((*itin1).coeff).
1044 add_dyn(ex_to_numeric((*itin2).coeff));
1045 if (expair_needs_further_processing(itin1))
1046 needs_further_processing = true;
1049 if (!ex_to_numeric((*itin1).coeff).is_zero()) {
1058 if (!ex_to_numeric((*itin1).coeff).is_zero()) {
1064 seq.erase(itout,last);
1067 if (needs_further_processing) {
1070 construct_from_epvector(v);
1075 #if EXPAIRSEQ_USE_HASHTAB
1077 unsigned expairseq::calc_hashtabsize(unsigned sz) const
1080 unsigned nearest_power_of_2 = 1 << log2(sz);
1081 // if (nearest_power_of_2 < maxhashtabsize/hashtabfactor) {
1082 // size = nearest_power_of_2*hashtabfactor;
1083 size = nearest_power_of_2/hashtabfactor;
1084 if (size<minhashtabsize)
1086 GINAC_ASSERT(hashtabsize<=0x8000000U); // really max size due to 31 bit hashing
1087 // hashtabsize must be a power of 2
1088 GINAC_ASSERT((1U << log2(size))==size);
1092 unsigned expairseq::calc_hashindex(const ex &e) const
1094 // calculate hashindex
1095 unsigned hash = e.gethash();
1097 if (is_a_numeric_hash(hash)) {
1098 hashindex = hashmask;
1100 hashindex = hash &hashmask;
1101 // last hashtab entry is reserved for numerics
1102 if (hashindex==hashmask) hashindex = 0;
1104 GINAC_ASSERT(hashindex>=0);
1105 GINAC_ASSERT((hashindex<hashtabsize)||(hashtabsize==0));
1109 void expairseq::shrink_hashtab(void)
1111 unsigned new_hashtabsize;
1112 while (hashtabsize!=(new_hashtabsize=calc_hashtabsize(seq.size()))) {
1113 GINAC_ASSERT(new_hashtabsize<hashtabsize);
1114 if (new_hashtabsize==0) {
1121 // shrink by a factor of 2
1122 unsigned half_hashtabsize = hashtabsize/2;
1123 for (unsigned i=0; i<half_hashtabsize-1; ++i)
1124 hashtab[i].merge(hashtab[i+half_hashtabsize],epp_is_less());
1125 // special treatment for numeric hashes
1126 hashtab[0].merge(hashtab[half_hashtabsize-1],epp_is_less());
1127 hashtab[half_hashtabsize-1] = hashtab[hashtabsize-1];
1128 hashtab.resize(half_hashtabsize);
1129 hashtabsize = half_hashtabsize;
1130 hashmask = hashtabsize-1;
1134 void expairseq::remove_hashtab_entry(epvector::const_iterator element)
1137 return; // nothing to do
1139 // calculate hashindex of element to be deleted
1140 unsigned hashindex = calc_hashindex((*element).rest);
1142 // find it in hashtab and remove it
1143 epplist &eppl = hashtab[hashindex];
1144 epplist::iterator epplit = eppl.begin();
1145 bool erased = false;
1146 while (epplit!=eppl.end()) {
1147 if (*epplit == element) {
1156 cout << "tried to erase " << element-seq.begin() << std::endl;
1157 cout << "size " << seq.end()-seq.begin() << std::endl;
1159 unsigned hashindex = calc_hashindex((*element).rest);
1160 epplist &eppl = hashtab[hashindex];
1161 epplist::iterator epplit=eppl.begin();
1163 while (epplit!=eppl.end()) {
1164 if (*epplit == element) {
1171 GINAC_ASSERT(erased);
1173 GINAC_ASSERT(erased);
1176 void expairseq::move_hashtab_entry(epvector::const_iterator oldpos,
1177 epvector::iterator newpos)
1179 GINAC_ASSERT(hashtabsize!=0);
1181 // calculate hashindex of element which was moved
1182 unsigned hashindex=calc_hashindex((*newpos).rest);
1184 // find it in hashtab and modify it
1185 epplist &eppl = hashtab[hashindex];
1186 epplist::iterator epplit = eppl.begin();
1187 while (epplit!=eppl.end()) {
1188 if (*epplit == oldpos) {
1194 GINAC_ASSERT(epplit!=eppl.end());
1197 void expairseq::sorted_insert(epplist &eppl, epp elem)
1199 epplist::iterator current = eppl.begin();
1200 while ((current!=eppl.end())&&((*(*current)).is_less(*elem))) {
1203 eppl.insert(current,elem);
1206 void expairseq::build_hashtab_and_combine(epvector::iterator &first_numeric,
1207 epvector::iterator &last_non_zero,
1208 vector<bool> &touched,
1209 unsigned &number_of_zeroes)
1211 epp current=seq.begin();
1213 while (current!=first_numeric) {
1214 if (is_ex_exactly_of_type((*current).rest,numeric)) {
1216 iter_swap(current,first_numeric);
1218 // calculate hashindex
1219 unsigned currenthashindex = calc_hashindex((*current).rest);
1221 // test if there is already a matching expair in the hashtab-list
1222 epplist &eppl=hashtab[currenthashindex];
1223 epplist::iterator epplit = eppl.begin();
1224 while (epplit!=eppl.end()) {
1225 if ((*current).rest.is_equal((*(*epplit)).rest))
1229 if (epplit==eppl.end()) {
1230 // no matching expair found, append this to end of list
1231 sorted_insert(eppl,current);
1234 // epplit points to a matching expair, combine it with current
1235 (*(*epplit)).coeff = ex_to_numeric((*(*epplit)).coeff).
1236 add_dyn(ex_to_numeric((*current).coeff));
1238 // move obsolete current expair to end by swapping with last_non_zero element
1239 // if this was a numeric, it is swapped with the expair before first_numeric
1240 iter_swap(current,last_non_zero);
1242 if (first_numeric!=last_non_zero) iter_swap(first_numeric,current);
1245 // test if combined term has coeff 0 and can be removed is done later
1246 touched[(*epplit)-seq.begin()]=true;
1252 void expairseq::drop_coeff_0_terms(epvector::iterator &first_numeric,
1253 epvector::iterator &last_non_zero,
1254 vector<bool> &touched,
1255 unsigned &number_of_zeroes)
1257 // move terms with coeff 0 to end and remove them from hashtab
1258 // check only those elements which have been touched
1259 epp current = seq.begin();
1261 while (current!=first_numeric) {
1265 } else if (!ex_to_numeric((*current).coeff).is_zero()) {
1269 remove_hashtab_entry(current);
1271 // move element to the end, unless it is already at the end
1272 if (current!=last_non_zero) {
1273 iter_swap(current,last_non_zero);
1275 bool numeric_swapped=first_numeric!=last_non_zero;
1276 if (numeric_swapped) iter_swap(first_numeric,current);
1277 epvector::iterator changed_entry;
1279 if (numeric_swapped)
1280 changed_entry = first_numeric;
1282 changed_entry = last_non_zero;
1287 if (first_numeric!=current) {
1289 // change entry in hashtab which referred to first_numeric or last_non_zero to current
1290 move_hashtab_entry(changed_entry,current);
1291 touched[current-seq.begin()] = touched[changed_entry-seq.begin()];
1300 GINAC_ASSERT(i==current-seq.begin());
1303 /** True if one of the coeffs vanishes, otherwise false.
1304 * This would be an invariant violation, so this should only be used for
1305 * debugging purposes. */
1306 bool expairseq::has_coeff_0(void) const
1308 for (epvector::const_iterator cit=seq.begin(); cit!=seq.end(); ++cit) {
1309 if ((*cit).coeff.is_zero())
1315 void expairseq::add_numerics_to_hashtab(epvector::iterator first_numeric,
1316 epvector::const_iterator last_non_zero)
1318 if (first_numeric==seq.end()) return; // no numerics
1320 epvector::iterator current = first_numeric;
1321 epvector::const_iterator last = last_non_zero+1;
1322 while (current!=last) {
1323 sorted_insert(hashtab[hashmask],current);
1328 void expairseq::combine_same_terms(void)
1330 // combine same terms, drop term with coeff 0, move numerics to end
1332 // calculate size of hashtab
1333 hashtabsize = calc_hashtabsize(seq.size());
1335 // hashtabsize is a power of 2
1336 hashmask = hashtabsize-1;
1340 hashtab.resize(hashtabsize);
1342 if (hashtabsize==0) {
1344 combine_same_terms_sorted_seq();
1345 GINAC_ASSERT(!has_coeff_0());
1349 // iterate through seq, move numerics to end,
1350 // fill hashtab and combine same terms
1351 epvector::iterator first_numeric = seq.end();
1352 epvector::iterator last_non_zero = seq.end()-1;
1354 vector<bool> touched;
1355 touched.reserve(seq.size());
1356 for (unsigned i=0; i<seq.size(); ++i) touched[i]=false;
1358 unsigned number_of_zeroes = 0;
1360 GINAC_ASSERT(!has_coeff_0());
1361 build_hashtab_and_combine(first_numeric,last_non_zero,touched,number_of_zeroes);
1363 cout << "in combine:" << std::endl;
1365 cout << "size=" << seq.end() - seq.begin() << std::endl;
1366 cout << "first_numeric=" << first_numeric - seq.begin() << std::endl;
1367 cout << "last_non_zero=" << last_non_zero - seq.begin() << std::endl;
1368 for (unsigned i=0; i<seq.size(); ++i) {
1369 if (touched[i]) cout << i << " is touched" << std::endl;
1371 cout << "end in combine" << std::endl;
1374 // there should not be any terms with coeff 0 from the beginning,
1375 // so it should be safe to skip this step
1376 if (number_of_zeroes!=0) {
1377 drop_coeff_0_terms(first_numeric,last_non_zero,touched,number_of_zeroes);
1379 cout << "in combine after drop:" << std::endl;
1381 cout << "size=" << seq.end() - seq.begin() << std::endl;
1382 cout << "first_numeric=" << first_numeric - seq.begin() << std::endl;
1383 cout << "last_non_zero=" << last_non_zero - seq.begin() << std::endl;
1384 for (unsigned i=0; i<seq.size(); ++i) {
1385 if (touched[i]) cout << i << " is touched" << std::endl;
1387 cout << "end in combine after drop" << std::endl;
1391 add_numerics_to_hashtab(first_numeric,last_non_zero);
1393 // pop zero elements
1394 for (unsigned i=0; i<number_of_zeroes; ++i) {
1398 // shrink hashtabsize to calculated value
1399 GINAC_ASSERT(!has_coeff_0());
1403 GINAC_ASSERT(!has_coeff_0());
1406 #endif // EXPAIRSEQ_USE_HASHTAB
1408 /** Check if this expairseq is in sorted (canonical) form. Useful mainly for
1409 * debugging or in assertions since being sorted is an invariance. */
1410 bool expairseq::is_canonical() const
1415 #if EXPAIRSEQ_USE_HASHTAB
1416 if (hashtabsize>0) return 1; // not canoncalized
1417 #endif // EXPAIRSEQ_USE_HASHTAB
1419 epvector::const_iterator it = seq.begin();
1420 epvector::const_iterator it_last = it;
1421 for (++it; it!=seq.end(); it_last=it, ++it) {
1422 if (!((*it_last).is_less(*it)||(*it_last).is_equal(*it))) {
1423 if (!is_ex_exactly_of_type((*it_last).rest,numeric)||
1424 !is_ex_exactly_of_type((*it).rest,numeric)) {
1425 // double test makes it easier to set a breakpoint...
1426 if (!is_ex_exactly_of_type((*it_last).rest,numeric)||
1427 !is_ex_exactly_of_type((*it).rest,numeric)) {
1428 printpair(std::clog,*it_last,0);
1430 printpair(std::clog,*it,0);
1432 std::clog << "pair1:" << std::endl;
1433 (*it_last).rest.printtree(std::clog);
1434 (*it_last).coeff.printtree(std::clog);
1435 std::clog << "pair2:" << std::endl;
1436 (*it).rest.printtree(std::clog);
1437 (*it).coeff.printtree(std::clog);
1447 /** Member-wise expand the expairs in this sequence.
1449 * @see expairseq::expand()
1450 * @return pointer to epvector containing expanded pairs or zero pointer,
1451 * if no members were changed. */
1452 epvector * expairseq::expandchildren(unsigned options) const
1454 epvector::const_iterator last = seq.end();
1455 epvector::const_iterator cit = seq.begin();
1457 const ex &expanded_ex = (*cit).rest.expand(options);
1458 if (!are_ex_trivially_equal((*cit).rest,expanded_ex)) {
1460 // something changed, copy seq, eval and return it
1461 epvector *s = new epvector;
1462 s->reserve(seq.size());
1464 // copy parts of seq which are known not to have changed
1465 epvector::const_iterator cit2 = seq.begin();
1467 s->push_back(*cit2);
1470 // copy first changed element
1471 s->push_back(combine_ex_with_coeff_to_pair(expanded_ex,
1475 while (cit2!=last) {
1476 s->push_back(combine_ex_with_coeff_to_pair((*cit2).rest.expand(options),
1485 return 0; // signalling nothing has changed
1489 /** Member-wise evaluate the expairs in this sequence.
1491 * @see expairseq::eval()
1492 * @return pointer to epvector containing evaluated pairs or zero pointer,
1493 * if no members were changed. */
1494 epvector * expairseq::evalchildren(int level) const
1496 // returns a NULL pointer if nothing had to be evaluated
1497 // returns a pointer to a newly created epvector otherwise
1498 // (which has to be deleted somewhere else)
1503 if (level == -max_recursion_level)
1504 throw(std::runtime_error("max recursion level reached"));
1507 epvector::const_iterator last=seq.end();
1508 epvector::const_iterator cit=seq.begin();
1510 const ex &evaled_ex = (*cit).rest.eval(level);
1511 if (!are_ex_trivially_equal((*cit).rest,evaled_ex)) {
1513 // something changed, copy seq, eval and return it
1514 epvector *s = new epvector;
1515 s->reserve(seq.size());
1517 // copy parts of seq which are known not to have changed
1518 epvector::const_iterator cit2=seq.begin();
1520 s->push_back(*cit2);
1523 // copy first changed element
1524 s->push_back(combine_ex_with_coeff_to_pair(evaled_ex,
1528 while (cit2!=last) {
1529 s->push_back(combine_ex_with_coeff_to_pair((*cit2).rest.eval(level),
1538 return 0; // signalling nothing has changed
1542 /** Member-wise evaluate numerically all expairs in this sequence.
1544 * @see expairseq::evalf()
1545 * @return epvector with all entries evaluated numerically. */
1546 epvector expairseq::evalfchildren(int level) const
1551 if (level==-max_recursion_level)
1552 throw(std::runtime_error("max recursion level reached"));
1555 s.reserve(seq.size());
1558 for (epvector::const_iterator it=seq.begin(); it!=seq.end(); ++it) {
1559 s.push_back(combine_ex_with_coeff_to_pair((*it).rest.evalf(level),
1560 (*it).coeff.evalf(level)));
1566 /** Member-wise normalize all expairs in this sequence.
1568 * @see expairseq::normal()
1569 * @return epvector with all entries normalized. */
1570 epvector expairseq::normalchildren(int level) const
1575 if (level==-max_recursion_level)
1576 throw(std::runtime_error("max recursion level reached"));
1579 s.reserve(seq.size());
1582 for (epvector::const_iterator it=seq.begin(); it!=seq.end(); ++it) {
1583 s.push_back(combine_ex_with_coeff_to_pair((*it).rest.normal(level),
1590 /** Member-wise differentiate all expairs in this sequence.
1592 * @see expairseq::diff()
1593 * @return epvector with all entries differentiated. */
1594 epvector expairseq::diffchildren(const symbol &y) const
1597 s.reserve(seq.size());
1599 for (epvector::const_iterator it=seq.begin(); it!=seq.end(); ++it) {
1600 s.push_back(combine_ex_with_coeff_to_pair((*it).rest.diff(y),
1607 /** Member-wise substitute in this sequence.
1609 * @see expairseq::subs()
1610 * @return pointer to epvector containing pairs after application of subs or zero
1611 * pointer, if no members were changed. */
1612 epvector * expairseq::subschildren(const lst &ls, const lst &lr) const
1614 // returns a NULL pointer if nothing had to be substituted
1615 // returns a pointer to a newly created epvector otherwise
1616 // (which has to be deleted somewhere else)
1617 GINAC_ASSERT(ls.nops()==lr.nops());
1619 epvector::const_iterator last = seq.end();
1620 epvector::const_iterator cit = seq.begin();
1622 const ex &subsed_ex=(*cit).rest.subs(ls,lr);
1623 if (!are_ex_trivially_equal((*cit).rest,subsed_ex)) {
1625 // something changed, copy seq, subs and return it
1626 epvector *s = new epvector;
1627 s->reserve(seq.size());
1629 // copy parts of seq which are known not to have changed
1630 epvector::const_iterator cit2 = seq.begin();
1632 s->push_back(*cit2);
1635 // copy first changed element
1636 s->push_back(combine_ex_with_coeff_to_pair(subsed_ex,
1640 while (cit2!=last) {
1641 s->push_back(combine_ex_with_coeff_to_pair((*cit2).rest.subs(ls,lr),
1650 return 0; // signalling nothing has changed
1654 // static member variables
1659 unsigned expairseq::precedence = 10;
1661 #if EXPAIRSEQ_USE_HASHTAB
1662 unsigned expairseq::maxhashtabsize = 0x4000000U;
1663 unsigned expairseq::minhashtabsize = 0x1000U;
1664 unsigned expairseq::hashtabfactor = 1;
1665 #endif // EXPAIRSEQ_USE_HASHTAB
1667 } // namespace GiNaC