]> www.ginac.de Git - ginac.git/blobdiff - ginac/expairseq.cpp
Fixed bug in ::match.
[ginac.git] / ginac / expairseq.cpp
index 24e2b5602d6e6dfe0c196f22ef220acdb3338133..9d72a7873c45ee23a4a4ceb0236a62571dfafe09 100644 (file)
@@ -3,7 +3,7 @@
  *  Implementation of sequences of expression pairs. */
 
 /*
- *  GiNaC Copyright (C) 1999-2003 Johannes Gutenberg University Mainz, Germany
+ *  GiNaC Copyright (C) 1999-2006 Johannes Gutenberg University Mainz, Germany
  *
  *  This program is free software; you can redistribute it and/or modify
  *  it under the terms of the GNU General Public License as published by
@@ -17,7 +17,7 @@
  *
  *  You should have received a copy of the GNU General Public License
  *  along with this program; if not, write to the Free Software
- *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+ *  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
  */
 
 #include <iostream>
@@ -27,6 +27,7 @@
 
 #include "expairseq.h"
 #include "lst.h"
+#include "add.h"
 #include "mul.h"
 #include "power.h"
 #include "relational.h"
@@ -34,6 +35,7 @@
 #include "archive.h"
 #include "operators.h"
 #include "utils.h"
+#include "indexed.h"
 
 #if EXPAIRSEQ_USE_HASHTAB
 #include <cmath>
@@ -66,7 +68,7 @@ public:
 
 // public
 
-expairseq::expairseq() : inherited(TINFO_expairseq)
+expairseq::expairseq() : inherited(&expairseq::tinfo_static)
 #if EXPAIRSEQ_USE_HASHTAB
                                                    , hashtabsize(0)
 #endif // EXPAIRSEQ_USE_HASHTAB
@@ -105,33 +107,32 @@ void expairseq::copy(const expairseq &other)
 // other constructors
 //////////
 
-expairseq::expairseq(const ex &lh, const ex &rh) : inherited(TINFO_expairseq)
+expairseq::expairseq(const ex &lh, const ex &rh) : inherited(&expairseq::tinfo_static)
 {
        construct_from_2_ex(lh,rh);
        GINAC_ASSERT(is_canonical());
 }
 
-expairseq::expairseq(const exvector &v) : inherited(TINFO_expairseq)
+expairseq::expairseq(const exvector &v) : inherited(&expairseq::tinfo_static)
 {
        construct_from_exvector(v);
        GINAC_ASSERT(is_canonical());
 }
 
 expairseq::expairseq(const epvector &v, const ex &oc)
-  : inherited(TINFO_expairseq), overall_coeff(oc)
+  : inherited(&expairseq::tinfo_static), overall_coeff(oc)
 {
        GINAC_ASSERT(is_a<numeric>(oc));
        construct_from_epvector(v);
        GINAC_ASSERT(is_canonical());
 }
 
-expairseq::expairseq(epvector *vp, const ex &oc)
-  : inherited(TINFO_expairseq), overall_coeff(oc)
+expairseq::expairseq(std::auto_ptr<epvector> vp, const ex &oc)
+  : inherited(&expairseq::tinfo_static), overall_coeff(oc)
 {
-       GINAC_ASSERT(vp!=0);
+       GINAC_ASSERT(vp.get()!=0);
        GINAC_ASSERT(is_a<numeric>(oc));
        construct_from_epvector(*vp);
-       delete vp;
        GINAC_ASSERT(is_canonical());
 }
 
@@ -283,7 +284,7 @@ ex expairseq::op(size_t i) const
 
 ex expairseq::map(map_function &f) const
 {
-       epvector *v = new epvector;
+       std::auto_ptr<epvector> v(new epvector);
        v->reserve(seq.size());
 
        epvector::const_iterator cit = seq.begin(), last = seq.end();
@@ -304,11 +305,58 @@ ex expairseq::eval(int level) const
        if ((level==1) && (flags &status_flags::evaluated))
                return *this;
        
-       epvector *vp = evalchildren(level);
-       if (vp==0)
+       std::auto_ptr<epvector> vp = evalchildren(level);
+       if (vp.get() == 0)
                return this->hold();
        
-       return (new expairseq(vp,overall_coeff))->setflag(status_flags::dynallocated | status_flags::evaluated);
+       return (new expairseq(vp, overall_coeff))->setflag(status_flags::dynallocated | status_flags::evaluated);
+}
+
+epvector* conjugateepvector(const epvector&epv)
+{
+       epvector *newepv = 0;
+       for (epvector::const_iterator i=epv.begin(); i!=epv.end(); ++i) {
+               if(newepv) {
+                       newepv->push_back(i->conjugate());
+                       continue;
+               }
+               expair x = i->conjugate();
+               if (x.is_equal(*i)) {
+                       continue;
+               }
+               newepv = new epvector;
+               newepv->reserve(epv.size());
+               for (epvector::const_iterator j=epv.begin(); j!=i; ++j) {
+                       newepv->push_back(*j);
+               }
+               newepv->push_back(x);
+       }
+       return newepv;
+}
+
+ex expairseq::conjugate() const
+{
+       epvector* newepv = conjugateepvector(seq);
+       ex x = overall_coeff.conjugate();
+       if (!newepv && are_ex_trivially_equal(x, overall_coeff)) {
+               return *this;
+       }
+       ex result = thisexpairseq(newepv ? *newepv : seq, x);
+       if (newepv) {
+               delete newepv;
+       }
+       return result;
+}
+
+bool expairseq::is_polynomial(const ex & var) const
+{
+       if (!is_exactly_a<add>(*this) && !is_exactly_a<mul>(*this))
+               return basic::is_polynomial(var);
+       for (epvector::const_iterator i=seq.begin(); i!=seq.end(); ++i) {
+               if (!(i->rest).is_polynomial(var))
+                       return false;
+       }
+       return true;
 }
 
 bool expairseq::match(const ex & pattern, lst & repl_lst) const
@@ -347,10 +395,20 @@ bool expairseq::match(const ex & pattern, lst & repl_lst) const
                                continue;
                        exvector::iterator it = ops.begin(), itend = ops.end();
                        while (it != itend) {
+                               lst::const_iterator last_el = repl_lst.end();
+                               --last_el;
                                if (it->match(p, repl_lst)) {
                                        ops.erase(it);
                                        goto found;
                                }
+                               while(true) {
+                                       lst::const_iterator next_el = last_el;
+                                       ++next_el;
+                                       if(next_el == repl_lst.end())
+                                               break;
+                                       else
+                                               repl_lst.remove_last();
+                               }
                                ++it;
                        }
                        return false; // no match found
@@ -363,7 +421,7 @@ found:              ;
                        // it has already been matched before, in which case the matches
                        // must be equal)
                        size_t num = ops.size();
-                       epvector *vp = new epvector();
+                       std::auto_ptr<epvector> vp(new epvector);
                        vp->reserve(num);
                        for (size_t i=0; i<num; i++)
                                vp->push_back(split_ex_to_pair(ops[i]));
@@ -387,8 +445,8 @@ found:              ;
 
 ex expairseq::subs(const exmap & m, unsigned options) const
 {
-       epvector *vp = subschildren(m, options);
-       if (vp)
+       std::auto_ptr<epvector> vp = subschildren(m, options);
+       if (vp.get())
                return ex_to<basic>(thisexpairseq(vp, overall_coeff));
        else if ((options & subs_options::algebraic) && is_exactly_a<mul>(*this))
                return static_cast<const mul *>(this)->algebraic_subs_mul(m, options);
@@ -535,7 +593,7 @@ unsigned expairseq::return_type() const
 
 unsigned expairseq::calchash() const
 {
-       unsigned v = golden_ratio_hash(this->tinfo());
+       unsigned v = golden_ratio_hash((p_int)this->tinfo());
        epvector::const_iterator i = seq.begin();
        const epvector::const_iterator end = seq.end();
        while (i != end) {
@@ -561,12 +619,13 @@ unsigned expairseq::calchash() const
 
 ex expairseq::expand(unsigned options) const
 {
-       epvector *vp = expandchildren(options);
-       if (vp == NULL) {
+       std::auto_ptr<epvector> vp = expandchildren(options);
+       if (vp.get())
+               return thisexpairseq(vp, overall_coeff);
+       else {
                // The terms have not changed, so it is safe to declare this expanded
                return (options == 0) ? setflag(status_flags::expanded) : *this;
-       } else
-               return thisexpairseq(vp, overall_coeff);
+       }
 }
 
 //////////
@@ -585,12 +644,12 @@ ex expairseq::expand(unsigned options) const
  *  definition. */
 ex expairseq::thisexpairseq(const epvector &v, const ex &oc) const
 {
-       return expairseq(v,oc);
+       return expairseq(v, oc);
 }
 
-ex expairseq::thisexpairseq(epvector *vp, const ex &oc) const
+ex expairseq::thisexpairseq(std::auto_ptr<epvector> vp, const ex &oc) const
 {
-       return expairseq(vp,oc);
+       return expairseq(vp, oc);
 }
 
 void expairseq::printpair(const print_context & c, const expair & p, unsigned upper_precedence) const
@@ -721,8 +780,15 @@ void expairseq::construct_from_2_ex(const ex &lh, const ex &rh)
                                construct_from_2_ex_via_exvector(lh,rh);
                        } else {
 #endif // EXPAIRSEQ_USE_HASHTAB
-                               construct_from_2_expairseq(ex_to<expairseq>(lh),
-                                                          ex_to<expairseq>(rh));
+                               if(is_a<mul>(lh))
+                               {       
+                                       ex newrh=rename_dummy_indices_uniquely(lh, rh);
+                                       construct_from_2_expairseq(ex_to<expairseq>(lh),
+                                                                  ex_to<expairseq>(newrh));
+                               }
+                               else
+                                       construct_from_2_expairseq(ex_to<expairseq>(lh),
+                                                                  ex_to<expairseq>(rh));
 #if EXPAIRSEQ_USE_HASHTAB
                        }
 #endif // EXPAIRSEQ_USE_HASHTAB
@@ -817,6 +883,7 @@ void expairseq::construct_from_2_expairseq(const expairseq &s1,
        
        while (first1!=last1 && first2!=last2) {
                int cmpval = (*first1).rest.compare((*first2).rest);
+
                if (cmpval==0) {
                        // combine terms
                        const numeric &newcoeff = ex_to<numeric>(first1->coeff).
@@ -972,13 +1039,27 @@ void expairseq::make_flat(const exvector &v)
        seq.reserve(v.size()+noperands-nexpairseqs);
        
        // copy elements and split off numerical part
+       exvector dummy_indices;
        cit = v.begin();
        while (cit!=v.end()) {
                if (ex_to<basic>(*cit).tinfo()==this->tinfo()) {
-                       const expairseq &subseqref = ex_to<expairseq>(*cit);
-                       combine_overall_coeff(subseqref.overall_coeff);
-                       epvector::const_iterator cit_s = subseqref.seq.begin();
-                       while (cit_s!=subseqref.seq.end()) {
+                       const expairseq *subseqref;
+                       ex newfactor;
+                       if(is_a<mul>(*cit))
+                       {
+                               exvector dummies_of_factor = get_all_dummy_indices(*cit);
+                               sort(dummies_of_factor.begin(), dummies_of_factor.end(), ex_is_less());
+                               newfactor = rename_dummy_indices_uniquely(dummy_indices, dummies_of_factor, *cit);
+                               subseqref = &(ex_to<expairseq>(newfactor));
+                               exvector new_dummy_indices;
+                               set_union(dummy_indices.begin(), dummy_indices.end(), dummies_of_factor.begin(), dummies_of_factor.end(), std::back_insert_iterator<exvector>(new_dummy_indices), ex_is_less());
+                               dummy_indices.swap(new_dummy_indices);
+                       }
+                       else
+                               subseqref = &ex_to<expairseq>(*cit);
+                       combine_overall_coeff(subseqref->overall_coeff);
+                       epvector::const_iterator cit_s = subseqref->seq.begin();
+                       while (cit_s!=subseqref->seq.end()) {
                                seq.push_back(*cit_s);
                                ++cit_s;
                        }
@@ -1450,7 +1531,7 @@ bool expairseq::is_canonical() const
  *  @see expairseq::expand()
  *  @return pointer to epvector containing expanded pairs or zero pointer,
  *  if no members were changed. */
-epvector * expairseq::expandchildren(unsigned options) const
+std::auto_ptr<epvector> expairseq::expandchildren(unsigned options) const
 {
        const epvector::const_iterator last = seq.end();
        epvector::const_iterator cit = seq.begin();
@@ -1459,7 +1540,7 @@ epvector * expairseq::expandchildren(unsigned options) const
                if (!are_ex_trivially_equal(cit->rest,expanded_ex)) {
                        
                        // something changed, copy seq, eval and return it
-                       epvector *s = new epvector;
+                       std::auto_ptr<epvector> s(new epvector);
                        s->reserve(seq.size());
                        
                        // copy parts of seq which are known not to have changed
@@ -1468,10 +1549,12 @@ epvector * expairseq::expandchildren(unsigned options) const
                                s->push_back(*cit2);
                                ++cit2;
                        }
+
                        // copy first changed element
                        s->push_back(combine_ex_with_coeff_to_pair(expanded_ex,
                                                                   cit2->coeff));
                        ++cit2;
+
                        // copy rest
                        while (cit2!=last) {
                                s->push_back(combine_ex_with_coeff_to_pair(cit2->rest.expand(options),
@@ -1483,7 +1566,7 @@ epvector * expairseq::expandchildren(unsigned options) const
                ++cit;
        }
        
-       return 0; // signalling nothing has changed
+       return std::auto_ptr<epvector>(0); // signalling nothing has changed
 }
 
 
@@ -1492,14 +1575,14 @@ epvector * expairseq::expandchildren(unsigned options) const
  *  @see expairseq::eval()
  *  @return pointer to epvector containing evaluated pairs or zero pointer,
  *  if no members were changed. */
-epvector * expairseq::evalchildren(int level) const
+std::auto_ptr<epvector> expairseq::evalchildren(int level) const
 {
        // returns a NULL pointer if nothing had to be evaluated
        // returns a pointer to a newly created epvector otherwise
        // (which has to be deleted somewhere else)
 
        if (level==1)
-               return 0;
+               return std::auto_ptr<epvector>(0);
        
        if (level == -max_recursion_level)
                throw(std::runtime_error("max recursion level reached"));
@@ -1512,7 +1595,7 @@ epvector * expairseq::evalchildren(int level) const
                if (!are_ex_trivially_equal(cit->rest,evaled_ex)) {
                        
                        // something changed, copy seq, eval and return it
-                       epvector *s = new epvector;
+                       std::auto_ptr<epvector> s(new epvector);
                        s->reserve(seq.size());
                        
                        // copy parts of seq which are known not to have changed
@@ -1521,10 +1604,12 @@ epvector * expairseq::evalchildren(int level) const
                                s->push_back(*cit2);
                                ++cit2;
                        }
+
                        // copy first changed element
                        s->push_back(combine_ex_with_coeff_to_pair(evaled_ex,
                                                                   cit2->coeff));
                        ++cit2;
+
                        // copy rest
                        while (cit2!=last) {
                                s->push_back(combine_ex_with_coeff_to_pair(cit2->rest.eval(level),
@@ -1536,16 +1621,77 @@ epvector * expairseq::evalchildren(int level) const
                ++cit;
        }
        
-       return 0; // signalling nothing has changed
+       return std::auto_ptr<epvector>(0); // signalling nothing has changed
+}
+
+class safe_inserter
+{
+       public:
+               safe_inserter(const ex&, const bool disable_renaming=false);
+               std::auto_ptr<epvector> getseq(){return epv;}
+               void insert_old_pair(const expair &p)
+               {
+                       epv->push_back(p);
+               }
+               void insert_new_pair(const expair &p, const ex &orig_ex);
+       private:
+               std::auto_ptr<epvector> epv;
+               bool dodummies;
+               exvector dummy_indices;
+               void update_dummy_indices(const exvector&);
+};
+
+safe_inserter::safe_inserter(const ex&e, const bool disable_renaming)
+               :epv(new epvector)
+{
+       epv->reserve(e.nops());
+       dodummies=is_a<mul>(e);
+       if(disable_renaming)
+               dodummies=false;
+       if(dodummies) {
+               dummy_indices = get_all_dummy_indices_safely(e);
+               sort(dummy_indices.begin(), dummy_indices.end(), ex_is_less());
+       }
+}
+
+void safe_inserter::update_dummy_indices(const exvector &v)
+{
+       exvector new_dummy_indices;
+       set_union(dummy_indices.begin(), dummy_indices.end(), v.begin(), v.end(),
+               std::back_insert_iterator<exvector>(new_dummy_indices), ex_is_less());
+       dummy_indices.swap(new_dummy_indices);
 }
 
+void safe_inserter::insert_new_pair(const expair &p, const ex &orig_ex)
+{
+       if(!dodummies) {
+               epv->push_back(p);
+               return;
+       }
+       exvector dummies_of_factor = get_all_dummy_indices_safely(p.rest);
+       if(dummies_of_factor.size() == 0) {
+               epv->push_back(p);
+               return;
+       }
+       sort(dummies_of_factor.begin(), dummies_of_factor.end(), ex_is_less());
+       exvector dummies_of_orig_ex = get_all_dummy_indices_safely(orig_ex);
+       sort(dummies_of_orig_ex.begin(), dummies_of_orig_ex.end(), ex_is_less());
+       exvector new_dummy_indices;
+       new_dummy_indices.reserve(dummy_indices.size());
+       set_difference(dummy_indices.begin(), dummy_indices.end(), dummies_of_orig_ex.begin(), dummies_of_orig_ex.end(),
+               std::back_insert_iterator<exvector>(new_dummy_indices), ex_is_less());
+       dummy_indices.swap(new_dummy_indices);
+       ex newfactor = rename_dummy_indices_uniquely(dummy_indices, dummies_of_factor, p.rest);
+       update_dummy_indices(dummies_of_factor);
+       epv -> push_back(expair(newfactor, p.coeff));
+}
 
 /** Member-wise substitute in this sequence.
  *
  *  @see expairseq::subs()
  *  @return pointer to epvector containing pairs after application of subs,
  *    or NULL pointer if no members were changed. */
-epvector * expairseq::subschildren(const exmap & m, unsigned options) const
+std::auto_ptr<epvector> expairseq::subschildren(const exmap & m, unsigned options) const
 {
        // When any of the objects to be substituted is a product or power
        // we have to recombine the pairs because the numeric coefficients may
@@ -1574,22 +1720,27 @@ epvector * expairseq::subschildren(const exmap & m, unsigned options) const
                        if (!are_ex_trivially_equal(orig_ex, subsed_ex)) {
 
                                // Something changed, copy seq, subs and return it
-                               epvector *s = new epvector;
-                               s->reserve(seq.size());
+                               safe_inserter s(*this, options & subs_options::no_index_renaming);
 
                                // Copy parts of seq which are known not to have changed
-                               s->insert(s->begin(), seq.begin(), cit);
+                               for(epvector::const_iterator i=seq.begin(); i!=cit; ++i)
+                                       s.insert_old_pair(*i);
 
                                // Copy first changed element
-                               s->push_back(split_ex_to_pair(subsed_ex));
+                               s.insert_new_pair(split_ex_to_pair(subsed_ex), orig_ex);
                                ++cit;
 
                                // Copy rest
                                while (cit != last) {
-                                       s->push_back(split_ex_to_pair(recombine_pair_to_ex(*cit).subs(m, options)));
+                                       ex orig_ex = recombine_pair_to_ex(*cit);
+                                       ex subsed_ex = orig_ex.subs(m, options);
+                                       if(are_ex_trivially_equal(orig_ex, subsed_ex))
+                                               s.insert_old_pair(*cit);
+                                       else
+                                               s.insert_new_pair(split_ex_to_pair(subsed_ex), orig_ex);
                                        ++cit;
                                }
-                               return s;
+                               return s.getseq();
                        }
 
                        ++cit;
@@ -1605,23 +1756,27 @@ epvector * expairseq::subschildren(const exmap & m, unsigned options) const
                        if (!are_ex_trivially_equal(cit->rest, subsed_ex)) {
                        
                                // Something changed, copy seq, subs and return it
-                               epvector *s = new epvector;
-                               s->reserve(seq.size());
+                               safe_inserter s(*this, options & subs_options::no_index_renaming);
 
                                // Copy parts of seq which are known not to have changed
-                               s->insert(s->begin(), seq.begin(), cit);
+                               for(epvector::const_iterator i=seq.begin(); i!=cit; ++i)
+                                       s.insert_old_pair(*i);
                        
                                // Copy first changed element
-                               s->push_back(combine_ex_with_coeff_to_pair(subsed_ex, cit->coeff));
+                               s.insert_new_pair(combine_ex_with_coeff_to_pair(subsed_ex, cit->coeff), cit->rest);
                                ++cit;
 
                                // Copy rest
                                while (cit != last) {
-                                       s->push_back(combine_ex_with_coeff_to_pair(cit->rest.subs(m, options),
-                                                                                  cit->coeff));
+                                       const ex &orig_ex = cit->rest;
+                                       const ex &subsed_ex = cit->rest.subs(m, options);
+                                       if(are_ex_trivially_equal(orig_ex, subsed_ex))
+                                               s.insert_old_pair(*cit);
+                                       else
+                                               s.insert_new_pair(combine_ex_with_coeff_to_pair(subsed_ex, cit->coeff), orig_ex);
                                        ++cit;
                                }
-                               return s;
+                               return s.getseq();
                        }
 
                        ++cit;
@@ -1629,7 +1784,7 @@ epvector * expairseq::subschildren(const exmap & m, unsigned options) const
        }
        
        // Nothing has changed
-       return NULL;
+       return std::auto_ptr<epvector>(0);
 }
 
 //////////