subs() and normal() use maps instead of lists, resulting in a huge performance
authorChristian Bauer <Christian.Bauer@uni-mainz.de>
Fri, 11 Jul 2003 19:28:43 +0000 (19:28 +0000)
committerChristian Bauer <Christian.Bauer@uni-mainz.de>
Fri, 11 Jul 2003 19:28:43 +0000 (19:28 +0000)
boost for subs()

26 files changed:
ginac/add.h
ginac/basic.cpp
ginac/basic.h
ginac/container.h
ginac/ex.cpp
ginac/ex.h
ginac/expairseq.cpp
ginac/expairseq.h
ginac/flags.h
ginac/idx.cpp
ginac/idx.h
ginac/matrix.cpp
ginac/matrix.h
ginac/mul.cpp
ginac/mul.h
ginac/normal.cpp
ginac/numeric.h
ginac/power.cpp
ginac/power.h
ginac/pseries.cpp
ginac/pseries.h
ginac/relational.cpp
ginac/relational.h
ginac/structure.h
ginac/symbol.h
ginac/tensor.cpp

index f2a94ed..20faea8 100644 (file)
@@ -54,7 +54,7 @@ public:
        ex eval(int level=0) const;
        ex evalm() const;
        ex series(const relational & r, int order, unsigned options = 0) const;
-       ex normal(lst &sym_lst, lst &repl_lst, int level=0) const;
+       ex normal(exmap & repl, int level=0) const;
        numeric integer_content() const;
        ex smod(const numeric &xi) const;
        numeric max_coefficient() const;
index 595d645..a296cf9 100644 (file)
@@ -524,22 +524,19 @@ bool basic::match(const ex & pattern, lst & repl_lst) const
 }
 
 /** Helper function for subs(). Does not recurse into subexpressions. */
-ex basic::subs_one_level(const lst & ls, const lst & lr, unsigned options) const
+ex basic::subs_one_level(const exmap & m, unsigned options) const
 {
-       GINAC_ASSERT(ls.nops() == lr.nops());
-
-       lst::const_iterator its, itr;
+       exmap::const_iterator it;
 
        if (options & subs_options::subs_no_pattern) {
-               for (its = ls.begin(), itr = lr.begin(); its != ls.end(); ++its, ++itr) {
-                       if (is_equal(ex_to<basic>(*its)))
-                               return *itr;
-               }
+               it = m.find(*this);
+               if (it != m.end())
+                       return it->second;
        } else {
-               for (its = ls.begin(), itr = lr.begin(); its != ls.end(); ++its, ++itr) {
+               for (it = m.begin(); it != m.end(); ++it) {
                        lst repl_lst;
-                       if (match(ex_to<basic>(*its), repl_lst))
-                               return itr->subs(repl_lst, options | subs_options::subs_no_pattern); // avoid infinite recursion when re-substituting the wildcards
+                       if (match(ex_to<basic>(it->first), repl_lst))
+                               return it->second.subs(repl_lst, options | subs_options::subs_no_pattern); // avoid infinite recursion when re-substituting the wildcards
                }
        }
 
@@ -548,7 +545,7 @@ ex basic::subs_one_level(const lst & ls, const lst & lr, unsigned options) const
 
 /** Substitute a set of objects by arbitrary expressions. The ex returned
  *  will already be evaluated. */
-ex basic::subs(const lst & ls, const lst & lr, unsigned options) const
+ex basic::subs(const exmap & m, unsigned options) const
 {
        size_t num = nops();
        if (num) {
@@ -556,7 +553,7 @@ ex basic::subs(const lst & ls, const lst & lr, unsigned options) const
                // Substitute in subexpressions
                for (size_t i=0; i<num; i++) {
                        const ex & orig_op = op(i);
-                       const ex & subsed_op = orig_op.subs(ls, lr, options);
+                       const ex & subsed_op = orig_op.subs(m, options);
                        if (!are_ex_trivially_equal(orig_op, subsed_op)) {
 
                                // Something changed, clone the object
@@ -569,16 +566,16 @@ ex basic::subs(const lst & ls, const lst & lr, unsigned options) const
 
                                // Substitute the other operands
                                for (; i<num; i++)
-                                       copy->let_op(i) = op(i).subs(ls, lr, options);
+                                       copy->let_op(i) = op(i).subs(m, options);
 
                                // Perform substitutions on the new object as a whole
-                               return copy->subs_one_level(ls, lr, options);
+                               return copy->subs_one_level(m, options);
                        }
                }
        }
 
        // Nothing changed or no subexpressions
-       return subs_one_level(ls, lr, options);
+       return subs_one_level(m, options);
 }
 
 /** Default interface of nth derivative ex::diff(s, n).  It should be called
@@ -737,35 +734,6 @@ ex basic::expand(unsigned options) const
 
 // public
 
-/** Substitute objects in an expression (syntactic substitution) and return
- *  the result as a new expression.  There are two valid types of
- *  replacement arguments: 1) a relational like object==ex and 2) a list of
- *  relationals lst(object1==ex1,object2==ex2,...), which is converted to
- *  subs(lst(object1,object2,...),lst(ex1,ex2,...)). */
-ex basic::subs(const ex & e, unsigned options) const
-{
-       if (e.info(info_flags::relation_equal)) {
-               return subs(lst(e.lhs()), lst(e.rhs()), options);
-       }
-       if (!e.info(info_flags::list)) {
-               throw(std::invalid_argument("basic::subs(ex): argument must be a list"));
-       }
-
-       // Split list into two
-       lst ls;
-       lst lr;
-       GINAC_ASSERT(is_a<lst>(e));
-       for (lst::const_iterator it = ex_to<lst>(e).begin(); it != ex_to<lst>(e).end(); ++it) {
-               ex r = *it;
-               if (!r.info(info_flags::relation_equal)) {
-                       throw(std::invalid_argument("basic::subs(ex): argument must be a list of equations"));
-               }
-               ls.append(r.op(0));
-               lr.append(r.op(1));
-       }
-       return subs(ls, lr, options);
-}
-
 /** Compare objects syntactically to establish canonical ordering.
  *  All compare functions return: -1 for *this less than other, 0 equal,
  *  1 greater. */
index a68c2ee..a5af323 100644 (file)
@@ -25,6 +25,7 @@
 
 #include <cstddef> // for size_t
 #include <vector>
+#include <map>
 // CINT needs <algorithm> to work properly with <vector>
 #include <algorithm>
 
@@ -45,6 +46,7 @@ class print_context;
 template <class> class ptr;
 
 typedef std::vector<ex> exvector;
+typedef std::map<ex, ex, ex_is_less> exmap;
 
 
 /** Function object for map(). */
@@ -135,7 +137,7 @@ protected:
 public:
 
        // substitutions
-       virtual ex subs(const lst & ls, const lst & lr, unsigned options = 0) const;
+       virtual ex subs(const exmap & m, unsigned options = 0) const;
 
        // function mapping
        virtual ex map(map_function & f) const;
@@ -163,7 +165,7 @@ public:
        virtual ex series(const relational & r, int order, unsigned options = 0) const;
 
        // rational functions
-       virtual ex normal(lst &sym_lst, lst &repl_lst, int level = 0) const;
+       virtual ex normal(exmap & repl, int level = 0) const;
        virtual ex to_rational(lst &repl_lst) const;
        virtual ex to_polynomial(lst &repl_lst) const;
 
@@ -190,8 +192,7 @@ protected: // functions that should be called from class ex only
        
        // non-virtual functions in this class
 public:
-       ex subs(const ex & e, unsigned options = 0) const;
-       ex subs_one_level(const lst & ls, const lst & lr, unsigned options) const;
+       ex subs_one_level(const exmap & m, unsigned options) const;
        ex diff(const symbol & s, unsigned nth = 1) const;
        int compare(const basic & other) const;
        bool is_equal(const basic & other) const;
index 135adb9..9a40f2e 100644 (file)
@@ -275,7 +275,7 @@ public:
        ex op(size_t i) const;
        ex & let_op(size_t i);
        ex eval(int level = 0) const;
-       ex subs(const lst & ls, const lst & lr, unsigned options = 0) const;
+       ex subs(const exmap & m, unsigned options = 0) const;
 
 protected:
        bool is_equal_same_type(const basic & other) const;
@@ -328,7 +328,7 @@ public:
 
 protected:
        STLT evalchildren(int level) const;
-       STLT *subschildren(const lst & ls, const lst & lr, unsigned options = 0) const;
+       STLT *subschildren(const exmap & m, unsigned options = 0) const;
 };
 
 /** Default constructor */
@@ -424,13 +424,13 @@ ex container<C>::eval(int level) const
 }
 
 template <template <class> class C>
-ex container<C>::subs(const lst & ls, const lst & lr, unsigned options) const
+ex container<C>::subs(const exmap & m, unsigned options) const
 {
-       STLT *vp = subschildren(ls, lr, options);
+       STLT *vp = subschildren(m, options);
        if (vp)
-               return ex_to<basic>(thiscontainer(vp)).subs_one_level(ls, lr, options);
+               return ex_to<basic>(thiscontainer(vp)).subs_one_level(m, options);
        else
-               return subs_one_level(ls, lr, options);
+               return subs_one_level(m, options);
 }
 
 /** Compare two containers of the same type. */
@@ -587,7 +587,7 @@ typename container<C>::STLT container<C>::evalchildren(int level) const
 }
 
 template <template <class> class C>
-typename container<C>::STLT *container<C>::subschildren(const lst & ls, const lst & lr, unsigned options) const
+typename container<C>::STLT *container<C>::subschildren(const exmap & m, unsigned options) const
 {
        // returns a NULL pointer if nothing had to be substituted
        // returns a pointer to a newly created epvector otherwise
@@ -595,7 +595,7 @@ typename container<C>::STLT *container<C>::subschildren(const lst & ls, const ls
 
        const_iterator cit = seq.begin(), end = seq.end();
        while (cit != end) {
-               const ex & subsed_ex = cit->subs(ls, lr, options);
+               const ex & subsed_ex = cit->subs(m, options);
                if (!are_ex_trivially_equal(*cit, subsed_ex)) {
 
                        // copy first part of seq which hasn't changed
@@ -608,7 +608,7 @@ typename container<C>::STLT *container<C>::subschildren(const lst & ls, const ls
 
                        // copy rest
                        while (cit != end) {
-                               s->push_back(cit->subs(ls, lr, options));
+                               s->push_back(cit->subs(m, options));
                                ++cit;
                        }
 
index 43b9376..e53cbd3 100644 (file)
@@ -117,6 +117,65 @@ bool ex::find(const ex & pattern, lst & found) const
        return any_found;
 }
 
+/** Substitute objects in an expression (syntactic substitution) and return
+ *  the result as a new expression. */
+ex ex::subs(const lst & ls, const lst & lr, unsigned options) const
+{
+       GINAC_ASSERT(ls.nops() == lr.nops());
+
+       // Convert the lists to a map
+       exmap m;
+       for (lst::const_iterator its = ls.begin(), itr = lr.begin(); its != ls.end(); ++its, ++itr) {
+               m[*its] = *itr;
+
+               // Search for products and powers in the expressions to be substituted
+               // (for an optimization in expairseq::subs())
+               if (is_exactly_a<mul>(*its) || is_exactly_a<power>(*its))
+                       options |= subs_options::pattern_is_product;
+       }
+       if (!(options & subs_options::pattern_is_product))
+               options |= subs_options::pattern_is_not_product;
+
+       return bp->subs(m, options);
+}
+
+/** Substitute objects in an expression (syntactic substitution) and return
+ *  the result as a new expression.  There are two valid types of
+ *  replacement arguments: 1) a relational like object==ex and 2) a list of
+ *  relationals lst(object1==ex1,object2==ex2,...). */
+ex ex::subs(const ex & e, unsigned options) const
+{
+       if (e.info(info_flags::relation_equal)) {
+               exmap m;
+               const ex & s = e.lhs();
+               m[s] = e.rhs();
+               if (is_exactly_a<mul>(s) || is_exactly_a<power>(s))
+                       options |= subs_options::pattern_is_product;
+               return bp->subs(m, options);
+       } else if (!e.info(info_flags::list))
+               throw(std::invalid_argument("basic::subs(ex): argument must be a list"));
+
+       // Convert the list to a map
+       exmap m;
+       GINAC_ASSERT(is_a<lst>(e));
+       for (lst::const_iterator it = ex_to<lst>(e).begin(); it != ex_to<lst>(e).end(); ++it) {
+               ex r = *it;
+               if (!r.info(info_flags::relation_equal))
+                       throw(std::invalid_argument("basic::subs(ex): argument must be a list of equations"));
+               const ex & s = r.lhs();
+               m[s] = r.rhs();
+
+               // Search for products and powers in the expressions to be substituted
+               // (for an optimization in expairseq::subs())
+               if (is_exactly_a<mul>(s) || is_exactly_a<power>(s))
+                       options |= subs_options::pattern_is_product;
+       }
+       if (!(options & subs_options::pattern_is_product))
+               options |= subs_options::pattern_is_not_product;
+
+       return bp->subs(m, options);
+}
+
 /** Traverse expression tree with given visitor, preorder traversal. */
 void ex::traverse_preorder(visitor & v) const
 {
index eaa6b80..b936041 100644 (file)
@@ -131,8 +131,9 @@ public:
        bool match(const ex & pattern, lst & repl_lst) const { return bp->match(pattern, repl_lst); }
 
        // substitutions
-       ex subs(const lst & ls, const lst & lr, unsigned options = 0) const { return bp->subs(ls, lr, options); }
-       ex subs(const ex & e, unsigned options = 0) const { return bp->subs(e, options); }
+       ex subs(const exmap & m, unsigned options = 0) const;
+       ex subs(const lst & ls, const lst & lr, unsigned options = 0) const;
+       ex subs(const ex & e, unsigned options = 0) const;
 
        // function mapping
        ex map(map_function & f) const { return bp->map(f); }
@@ -429,12 +430,6 @@ inline ex series(const ex & thisex, const ex & r, int order, unsigned options =
 inline bool match(const ex & thisex, const ex & pattern, lst & repl_lst)
 { return thisex.match(pattern, repl_lst); }
 
-inline ex subs(const ex & thisex, const ex & e, unsigned options = 0)
-{ return thisex.subs(e, options); }
-
-inline ex subs(const ex & thisex, const lst & ls, const lst & lr, unsigned options = 0)
-{ return thisex.subs(ls, lr, options); }
-
 inline ex simplify_indexed(const ex & thisex)
 { return thisex.simplify_indexed(); }
 
@@ -497,13 +492,27 @@ struct ex_swap : public std::binary_function<ex, ex, void> {
        void operator() (ex &lh, ex &rh) const { lh.swap(rh); }
 };
 
+inline ex ex::subs(const exmap & m, unsigned options) const
+{
+       return bp->subs(m, options);
+}
+
+inline ex subs(const ex & thisex, const exmap & m, unsigned options = 0)
+{ return thisex.subs(m, options); }
+
+inline ex subs(const ex & thisex, const lst & ls, const lst & lr, unsigned options = 0)
+{ return thisex.subs(ls, lr, options); }
+
+inline ex subs(const ex & thisex, const ex & e, unsigned options = 0)
+{ return thisex.subs(e, options); }
+
 
 /* Convert function pointer to function object suitable for map(). */
 class pointer_to_map_function : public map_function {
 protected:
        ex (*ptr)(const ex &);
 public:
-       explicit pointer_to_map_function(ex (*x)(const ex &)) : ptr(x) {}
+       explicit pointer_to_map_function(ex x(const ex &)) : ptr(x) {}
        ex operator()(const ex & e) { return ptr(e); }
 };
 
@@ -513,7 +522,7 @@ protected:
        ex (*ptr)(const ex &, T1);
        T1 arg1;
 public:
-       explicit pointer_to_map_function_1arg(ex (*x)(const ex &, T1), T1 a1) : ptr(x), arg1(a1) {}
+       explicit pointer_to_map_function_1arg(ex x(const ex &, T1), T1 a1) : ptr(x), arg1(a1) {}
        ex operator()(const ex & e) { return ptr(e, arg1); }
 };
 
@@ -524,7 +533,7 @@ protected:
        T1 arg1;
        T2 arg2;
 public:
-       explicit pointer_to_map_function_2args(ex (*x)(const ex &, T1, T2), T1 a1, T2 a2) : ptr(x), arg1(a1), arg2(a2) {}
+       explicit pointer_to_map_function_2args(ex x(const ex &, T1, T2), T1 a1, T2 a2) : ptr(x), arg1(a1), arg2(a2) {}
        ex operator()(const ex & e) { return ptr(e, arg1, arg2); }
 };
 
@@ -536,11 +545,11 @@ protected:
        T2 arg2;
        T3 arg3;
 public:
-       explicit pointer_to_map_function_3args(ex (*x)(const ex &, T1, T2, T3), T1 a1, T2 a2, T3 a3) : ptr(x), arg1(a1), arg2(a2), arg3(a3) {}
+       explicit pointer_to_map_function_3args(ex x(const ex &, T1, T2, T3), T1 a1, T2 a2, T3 a3) : ptr(x), arg1(a1), arg2(a2), arg3(a3) {}
        ex operator()(const ex & e) { return ptr(e, arg1, arg2, arg3); }
 };
 
-inline ex ex::map(ex (*f)(const ex & e)) const
+inline ex ex::map(ex f(const ex &)) const
 {
        pointer_to_map_function fcn(f);
        return bp->map(fcn);
index f8675b7..5abdbc0 100644 (file)
@@ -386,15 +386,15 @@ found:            ;
        return inherited::match(pattern, repl_lst);
 }
 
-ex expairseq::subs(const lst &ls, const lst &lr, unsigned options) const
+ex expairseq::subs(const exmap & m, unsigned options) const
 {
-       epvector *vp = subschildren(ls, lr, options);
+       epvector *vp = subschildren(m, options);
        if (vp)
                return ex_to<basic>(thisexpairseq(vp, overall_coeff));
        else if ((options & subs_options::subs_algebraic) && is_exactly_a<mul>(*this))
-               return static_cast<const mul *>(this)->algebraic_subs_mul(ls, lr, options);
+               return static_cast<const mul *>(this)->algebraic_subs_mul(m, options);
        else
-               return subs_one_level(ls, lr, options);
+               return subs_one_level(m, options);
 }
 
 // protected
@@ -1546,29 +1546,32 @@ epvector * expairseq::evalchildren(int level) const
  *  @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 lst &ls, const lst &lr, unsigned options) const
+epvector * expairseq::subschildren(const exmap & m, unsigned options) const
 {
-       GINAC_ASSERT(ls.nops()==lr.nops());
-
-       // The substitution is "complex" when any of the objects to be substituted
-       // is a product or power. In this case we have to recombine the pairs
-       // because the numeric coefficients may be part of the search pattern.
-       bool complex_subs = false;
-       for (lst::const_iterator it = ls.begin(); it != ls.end(); ++it) {
-               if (is_exactly_a<mul>(*it) || is_exactly_a<power>(*it)) {
-                       complex_subs = true;
-                       break;
+       // When any of the objects to be substituted is a product or power
+       // we have to recombine the pairs because the numeric coefficients may
+       // be part of the search pattern.
+       if (!(options & (subs_options::pattern_is_product | subs_options::pattern_is_not_product))) {
+
+               // Search the list of substitutions and cache our findings
+               for (exmap::const_iterator it = m.begin(); it != m.end(); ++it) {
+                       if (is_exactly_a<mul>(it->first) || is_exactly_a<power>(it->first)) {
+                               options |= subs_options::pattern_is_product;
+                               break;
+                       }
                }
+               if (!(options & subs_options::pattern_is_product))
+                       options |= subs_options::pattern_is_not_product;
        }
 
-       if (complex_subs) {
+       if (options & subs_options::pattern_is_product) {
 
                // Substitute in the recombined pairs
                epvector::const_iterator cit = seq.begin(), last = seq.end();
                while (cit != last) {
 
                        const ex &orig_ex = recombine_pair_to_ex(*cit);
-                       const ex &subsed_ex = orig_ex.subs(ls, lr, options);
+                       const ex &subsed_ex = orig_ex.subs(m, options);
                        if (!are_ex_trivially_equal(orig_ex, subsed_ex)) {
 
                                // Something changed, copy seq, subs and return it
@@ -1584,7 +1587,7 @@ epvector * expairseq::subschildren(const lst &ls, const lst &lr, unsigned option
 
                                // Copy rest
                                while (cit != last) {
-                                       s->push_back(split_ex_to_pair(recombine_pair_to_ex(*cit).subs(ls, lr, options)));
+                                       s->push_back(split_ex_to_pair(recombine_pair_to_ex(*cit).subs(m, options)));
                                        ++cit;
                                }
                                return s;
@@ -1599,7 +1602,7 @@ epvector * expairseq::subschildren(const lst &ls, const lst &lr, unsigned option
                epvector::const_iterator cit = seq.begin(), last = seq.end();
                while (cit != last) {
 
-                       const ex &subsed_ex = cit->rest.subs(ls, lr, options);
+                       const ex &subsed_ex = cit->rest.subs(m, options);
                        if (!are_ex_trivially_equal(cit->rest, subsed_ex)) {
                        
                                // Something changed, copy seq, subs and return it
@@ -1615,7 +1618,7 @@ epvector * expairseq::subschildren(const lst &ls, const lst &lr, unsigned option
 
                                // Copy rest
                                while (cit != last) {
-                                       s->push_back(combine_ex_with_coeff_to_pair(cit->rest.subs(ls, lr, options),
+                                       s->push_back(combine_ex_with_coeff_to_pair(cit->rest.subs(m, options),
                                                                                   cit->coeff));
                                        ++cit;
                                }
index f108635..b6f9bc7 100644 (file)
@@ -80,7 +80,7 @@ public:
        ex to_rational(lst &repl_lst) const;
        ex to_polynomial(lst &repl_lst) const;
        bool match(const ex & pattern, lst & repl_lst) const;
-       ex subs(const lst & ls, const lst & lr, unsigned options = 0) const;
+       ex subs(const exmap & m, unsigned options = 0) const;
 protected:
        bool is_equal_same_type(const basic & other) const;
        unsigned return_type() const;
@@ -146,7 +146,7 @@ protected:
        bool is_canonical() const;
        epvector * expandchildren(unsigned options) const;
        epvector * evalchildren(int level) const;
-       epvector * subschildren(const lst & ls, const lst & lr, unsigned options = 0) const;
+       epvector * subschildren(const exmap & m, unsigned options = 0) const;
        
 // member variables
        
index 04c430c..c14e584 100644 (file)
@@ -29,9 +29,8 @@ namespace GiNaC {
 class expand_options {
 public:
        enum {
-               expand_trigonometric = 0x0001,
-               expand_indexed = 0x0002,
-               expand_function_args = 0x0004
+               expand_indexed = 0x0001,      ///< expands (a+b).i to a.i+b.i
+               expand_function_args = 0x0002 ///< expands the arguments of functions
        };
 };
 
@@ -39,8 +38,12 @@ public:
 class subs_options {
 public:
        enum {
-               subs_no_pattern = 0x0001,
-               subs_algebraic = 0x0002
+               // this should be called "no_pattern"...
+               subs_no_pattern = 0x0001,        ///< disable pattern matching
+               // this should be called "algebraic"...
+               subs_algebraic = 0x0002,         ///< enable algebraic substitutions
+        pattern_is_product = 0x0004,     ///< used internally by expairseq::subschildren()
+        pattern_is_not_product = 0x0008  ///< used internally by expairseq::subschildren()
        };
 };
 
@@ -163,10 +166,10 @@ public:
 class status_flags {
 public:
        enum {
-               dynallocated    = 0x0001,       ///< Heap-allocated (i.e. created by new if we want to be clever and bypass the stack, @see ex::construct_from_basic() )
-               evaluated       = 0x0002,       ///< .eval() has already done its job
-               expanded        = 0x0004,       ///< .expand(0) has already done its job (other expand() options ignore this flag)
-               hash_calculated = 0x0008        ///< .calchash() has already done its job
+               dynallocated    = 0x0001, ///< heap-allocated (i.e. created by new if we want to be clever and bypass the stack, @see ex::construct_from_basic() )
+               evaluated       = 0x0002, ///< .eval() has already done its job
+               expanded        = 0x0004, ///< .expand(0) has already done its job (other expand() options ignore this flag)
+               hash_calculated = 0x0008  ///< .calchash() has already done its job
        };
 };
 
index 265c107..cc5f52c 100644 (file)
@@ -348,29 +348,25 @@ ex idx::evalf(int level) const
        return *this;
 }
 
-ex idx::subs(const lst & ls, const lst & lr, unsigned options) const
+ex idx::subs(const exmap & m, unsigned options) const
 {
-       GINAC_ASSERT(ls.nops() == lr.nops());
-
        // First look for index substitutions
-       lst::const_iterator its, itr;
-       for (its = ls.begin(), itr = lr.begin(); its != ls.end(); ++its, ++itr) {
-               if (is_equal(ex_to<basic>(*its))) {
-
-                       // Substitution index->index
-                       if (is_a<idx>(*itr))
-                               return *itr;
-
-                       // Otherwise substitute value
-                       idx *i_copy = duplicate();
-                       i_copy->value = *itr;
-                       i_copy->clearflag(status_flags::hash_calculated);
-                       return i_copy->setflag(status_flags::dynallocated);
-               }
+       exmap::const_iterator it = m.find(*this);
+       if (it != m.end()) {
+
+               // Substitution index->index
+               if (is_a<idx>(it->second))
+                       return it->second;
+
+               // Otherwise substitute value
+               idx *i_copy = duplicate();
+               i_copy->value = it->second;
+               i_copy->clearflag(status_flags::hash_calculated);
+               return i_copy->setflag(status_flags::dynallocated);
        }
 
        // None, substitute objects in value (not in dimension)
-       const ex &subsed_value = value.subs(ls, lr, options);
+       const ex &subsed_value = value.subs(m, options);
        if (are_ex_trivially_equal(value, subsed_value))
                return *this;
 
index bb74dda..b6e63e0 100644 (file)
@@ -54,7 +54,7 @@ public:
        ex op(size_t i) const;
        ex map(map_function & f) const;
        ex evalf(int level = 0) const;
-       ex subs(const lst & ls, const lst & lr, unsigned options = 0) const;
+       ex subs(const exmap & m, unsigned options = 0) const;
 
 protected:
        ex derivative(const symbol & s) const;
index 990da83..cdb875d 100644 (file)
@@ -229,14 +229,14 @@ ex matrix::eval(int level) const
                                                                                           status_flags::evaluated);
 }
 
-ex matrix::subs(const lst & ls, const lst & lr, unsigned options) const
-{      
+ex matrix::subs(const exmap & mp, unsigned options) const
+{
        exvector m2(row * col);
        for (unsigned r=0; r<row; ++r)
                for (unsigned c=0; c<col; ++c)
-                       m2[r*col+c] = m[r*col+c].subs(ls, lr, options);
+                       m2[r*col+c] = m[r*col+c].subs(mp, options);
 
-       return matrix(row, col, m2).subs_one_level(ls, lr, options);
+       return matrix(row, col, m2).subs_one_level(mp, options);
 }
 
 // protected
index d0d9d7d..f951850 100644 (file)
@@ -49,7 +49,7 @@ public:
        ex & let_op(size_t i);
        ex eval(int level=0) const;
        ex evalm() const {return *this;}
-       ex subs(const lst & ls, const lst & lr, unsigned options = 0) const;
+       ex subs(const exmap & m, unsigned options = 0) const;
        ex eval_indexed(const basic & i) const;
        ex add_indexed(const ex & self, const ex & other) const;
        ex scalar_mul_indexed(const ex & self, const numeric & other) const;
index 8fe7399..2c84565 100644 (file)
@@ -558,27 +558,26 @@ bool tryfactsubs(const ex & origfactor, const ex & patternfactor, int & nummatch
        return true;
 }
 
-ex mul::algebraic_subs_mul(const lst & ls, const lst & lr, unsigned options) const
+ex mul::algebraic_subs_mul(const exmap & m, unsigned options) const
 {      
        std::vector<bool> subsed(seq.size(), false);
        exvector subsresult(seq.size());
 
-       lst::const_iterator its, itr;
-       for (its = ls.begin(), itr = lr.begin(); its != ls.end(); ++its, ++itr) {
+       for (exmap::const_iterator it = m.begin(); it != m.end(); ++it) {
 
-               if (is_exactly_a<mul>(*its)) {
+               if (is_exactly_a<mul>(it->first)) {
 
                        int nummatches = std::numeric_limits<int>::max();
                        std::vector<bool> currsubsed(seq.size(), false);
                        bool succeed = true;
                        lst repls;
 
-                       for (size_t j=0; j<its->nops(); j++) {
+                       for (size_t j=0; j<it->first.nops(); j++) {
                                bool found=false;
                                for (size_t k=0; k<nops(); k++) {
                                        if (currsubsed[k] || subsed[k])
                                                continue;
-                                       if (tryfactsubs(op(k), its->op(j), nummatches, repls)) {
+                                       if (tryfactsubs(op(k), it->first.op(j), nummatches, repls)) {
                                                currsubsed[k] = true;
                                                found = true;
                                                break;
@@ -599,7 +598,7 @@ ex mul::algebraic_subs_mul(const lst & ls, const lst & lr, unsigned options) con
                                                subsresult[j] = op(j);
                                        else {
                                                foundfirstsubsedfactor = true;
-                                               subsresult[j] = op(j) * power(itr->subs(ex(repls), subs_options::subs_no_pattern) / its->subs(ex(repls), subs_options::subs_no_pattern), nummatches);
+                                               subsresult[j] = op(j) * power(it->second.subs(ex(repls), subs_options::subs_no_pattern) / it->first.subs(ex(repls), subs_options::subs_no_pattern), nummatches);
                                        }
                                        subsed[j] = true;
                                }
@@ -611,9 +610,9 @@ ex mul::algebraic_subs_mul(const lst & ls, const lst & lr, unsigned options) con
                        lst repls;
 
                        for (size_t j=0; j<this->nops(); j++) {
-                               if (!subsed[j] && tryfactsubs(op(j), *its, nummatches, repls)) {
+                               if (!subsed[j] && tryfactsubs(op(j), it->first, nummatches, repls)) {
                                        subsed[j] = true;
-                                       subsresult[j] = op(j) * power(itr->subs(ex(repls), subs_options::subs_no_pattern) / its->subs(ex(repls), subs_options::subs_no_pattern), nummatches);
+                                       subsresult[j] = op(j) * power(it->second.subs(ex(repls), subs_options::subs_no_pattern) / it->first.subs(ex(repls), subs_options::subs_no_pattern), nummatches);
                                }
                        }
                }
@@ -627,7 +626,7 @@ ex mul::algebraic_subs_mul(const lst & ls, const lst & lr, unsigned options) con
                }
        }
        if (!subsfound)
-               return subs_one_level(ls, lr, options | subs_options::subs_algebraic);
+               return subs_one_level(m, options | subs_options::subs_algebraic);
 
        exvector ev; ev.reserve(nops());
        for (size_t i=0; i<nops(); i++) {
index 01fb205..0ec5720 100644 (file)
@@ -57,7 +57,7 @@ public:
        ex evalf(int level=0) const;
        ex evalm() const;
        ex series(const relational & s, int order, unsigned options = 0) const;
-       ex normal(lst &sym_lst, lst &repl_lst, int level = 0) const;
+       ex normal(exmap & repl, int level = 0) const;
        numeric integer_content() const;
        ex smod(const numeric &xi) const;
        numeric max_coefficient() const;
@@ -85,7 +85,7 @@ protected:
        
        // non-virtual functions in this class
 public:
-       ex algebraic_subs_mul(const lst & ls, const lst & lr, unsigned options) const;
+       ex algebraic_subs_mul(const exmap & m, unsigned options) const;
 protected:
        epvector * expandchildren(unsigned options) const;
 };
index cd2ad8e..3a8a82f 100644 (file)
@@ -1681,25 +1681,22 @@ ex sqrfree_parfrac(const ex & a, const symbol & x)
 
 
 /** Create a symbol for replacing the expression "e" (or return a previously
- *  assigned symbol). The symbol is appended to sym_lst and returned, the
- *  expression is appended to repl_lst.
+ *  assigned symbol). The symbol and expression are appended to repl, for
+ *  a later application of subs().
  *  @see ex::normal */
-static ex replace_with_symbol(const ex &e, lst &sym_lst, lst &repl_lst)
+static ex replace_with_symbol(const ex & e, exmap & repl)
 {
-       // Expression already in repl_lst? Then return the assigned symbol
-       lst::const_iterator its, itr;
-       for (its = sym_lst.begin(), itr = repl_lst.begin(); itr != repl_lst.end(); ++its, ++itr)
-               if (itr->is_equal(e))
-                       return *its;
+       // Expression already in repl? Then return the assigned symbol
+       for (exmap::const_iterator it = repl.begin(); it != repl.end(); ++it)
+               if (it->second.is_equal(e))
+                       return it->first;
        
        // Otherwise create new symbol and add to list, taking care that the
-       // replacement expression doesn't contain symbols from the sym_lst
+       // replacement expression doesn't itself contain symbols from repl,
        // because subs() is not recursive
-       symbol s;
-       ex es(s);
-       ex e_replaced = e.subs(sym_lst, repl_lst);
-       sym_lst.append(es);
-       repl_lst.append(e_replaced);
+       ex es = (new symbol)->setflag(status_flags::dynallocated);
+       ex e_replaced = e.subs(repl);
+       repl[es] = e_replaced;
        return es;
 }
 
@@ -1708,7 +1705,7 @@ static ex replace_with_symbol(const ex &e, lst &sym_lst, lst &repl_lst)
  *  to repl_lst and the symbol is returned.
  *  @see basic::to_rational
  *  @see basic::to_polynomial */
-static ex replace_with_symbol(const ex &e, lst &repl_lst)
+static ex replace_with_symbol(const ex & e, lst & repl_lst)
 {
        // Expression already in repl_lst? Then return the assigned symbol
        for (lst::const_iterator it = repl_lst.begin(); it != repl_lst.end(); ++it)
@@ -1716,10 +1713,9 @@ static ex replace_with_symbol(const ex &e, lst &repl_lst)
                        return it->op(0);
        
        // Otherwise create new symbol and add to list, taking care that the
-       // replacement expression doesn't contain symbols from the sym_lst
+       // replacement expression doesn't itself contain symbols from the repl_lst,
        // because subs() is not recursive
-       symbol s;
-       ex es(s);
+       ex es = (new symbol)->setflag(status_flags::dynallocated);
        ex e_replaced = e.subs(repl_lst);
        repl_lst.append(es == e_replaced);
        return es;
@@ -1736,18 +1732,18 @@ struct normal_map_function : public map_function {
 /** Default implementation of ex::normal(). It normalizes the children and
  *  replaces the object with a temporary symbol.
  *  @see ex::normal */
-ex basic::normal(lst &sym_lst, lst &repl_lst, int level) const
+ex basic::normal(exmap & repl, int level) const
 {
        if (nops() == 0)
-               return (new lst(replace_with_symbol(*this, sym_lst, repl_lst), _ex1))->setflag(status_flags::dynallocated);
+               return (new lst(replace_with_symbol(*this, repl), _ex1))->setflag(status_flags::dynallocated);
        else {
                if (level == 1)
-                       return (new lst(replace_with_symbol(*this, sym_lst, repl_lst), _ex1))->setflag(status_flags::dynallocated);
+                       return (new lst(replace_with_symbol(*this, repl), _ex1))->setflag(status_flags::dynallocated);
                else if (level == -max_recursion_level)
                        throw(std::runtime_error("max recursion level reached"));
                else {
                        normal_map_function map_normal(level - 1);
-                       return (new lst(replace_with_symbol(map(map_normal), sym_lst, repl_lst), _ex1))->setflag(status_flags::dynallocated);
+                       return (new lst(replace_with_symbol(map(map_normal), repl), _ex1))->setflag(status_flags::dynallocated);
                }
        }
 }
@@ -1755,7 +1751,7 @@ ex basic::normal(lst &sym_lst, lst &repl_lst, int level) const
 
 /** Implementation of ex::normal() for symbols. This returns the unmodified symbol.
  *  @see ex::normal */
-ex symbol::normal(lst &sym_lst, lst &repl_lst, int level) const
+ex symbol::normal(exmap & repl, int level) const
 {
        return (new lst(*this, _ex1))->setflag(status_flags::dynallocated);
 }
@@ -1765,19 +1761,19 @@ ex symbol::normal(lst &sym_lst, lst &repl_lst, int level) const
  *  into re+I*im and replaces I and non-rational real numbers with a temporary
  *  symbol.
  *  @see ex::normal */
-ex numeric::normal(lst &sym_lst, lst &repl_lst, int level) const
+ex numeric::normal(exmap & repl, int level) const
 {
        numeric num = numer();
        ex numex = num;
 
        if (num.is_real()) {
                if (!num.is_integer())
-                       numex = replace_with_symbol(numex, sym_lst, repl_lst);
+                       numex = replace_with_symbol(numex, repl);
        } else { // complex
                numeric re = num.real(), im = num.imag();
-               ex re_ex = re.is_rational() ? re : replace_with_symbol(re, sym_lst, repl_lst);
-               ex im_ex = im.is_rational() ? im : replace_with_symbol(im, sym_lst, repl_lst);
-               numex = re_ex + im_ex * replace_with_symbol(I, sym_lst, repl_lst);
+               ex re_ex = re.is_rational() ? re : replace_with_symbol(re, repl);
+               ex im_ex = im.is_rational() ? im : replace_with_symbol(im, repl);
+               numex = re_ex + im_ex * replace_with_symbol(I, repl);
        }
 
        // Denominator is always a real integer (see numeric::denom())
@@ -1849,10 +1845,10 @@ static ex frac_cancel(const ex &n, const ex &d)
 /** Implementation of ex::normal() for a sum. It expands terms and performs
  *  fractional addition.
  *  @see ex::normal */
-ex add::normal(lst &sym_lst, lst &repl_lst, int level) const
+ex add::normal(exmap & repl, int level) const
 {
        if (level == 1)
-               return (new lst(replace_with_symbol(*this, sym_lst, repl_lst), _ex1))->setflag(status_flags::dynallocated);
+               return (new lst(replace_with_symbol(*this, repl), _ex1))->setflag(status_flags::dynallocated);
        else if (level == -max_recursion_level)
                throw(std::runtime_error("max recursion level reached"));
 
@@ -1862,12 +1858,12 @@ ex add::normal(lst &sym_lst, lst &repl_lst, int level) const
        dens.reserve(seq.size()+1);
        epvector::const_iterator it = seq.begin(), itend = seq.end();
        while (it != itend) {
-               ex n = ex_to<basic>(recombine_pair_to_ex(*it)).normal(sym_lst, repl_lst, level-1);
+               ex n = ex_to<basic>(recombine_pair_to_ex(*it)).normal(repl, level-1);
                nums.push_back(n.op(0));
                dens.push_back(n.op(1));
                it++;
        }
-       ex n = ex_to<numeric>(overall_coeff).normal(sym_lst, repl_lst, level-1);
+       ex n = ex_to<numeric>(overall_coeff).normal(repl, level-1);
        nums.push_back(n.op(0));
        dens.push_back(n.op(1));
        GINAC_ASSERT(nums.size() == dens.size());
@@ -1908,10 +1904,10 @@ ex add::normal(lst &sym_lst, lst &repl_lst, int level) const
 /** Implementation of ex::normal() for a product. It cancels common factors
  *  from fractions.
  *  @see ex::normal() */
-ex mul::normal(lst &sym_lst, lst &repl_lst, int level) const
+ex mul::normal(exmap & repl, int level) const
 {
        if (level == 1)
-               return (new lst(replace_with_symbol(*this, sym_lst, repl_lst), _ex1))->setflag(status_flags::dynallocated);
+               return (new lst(replace_with_symbol(*this, repl), _ex1))->setflag(status_flags::dynallocated);
        else if (level == -max_recursion_level)
                throw(std::runtime_error("max recursion level reached"));
 
@@ -1921,12 +1917,12 @@ ex mul::normal(lst &sym_lst, lst &repl_lst, int level) const
        ex n;
        epvector::const_iterator it = seq.begin(), itend = seq.end();
        while (it != itend) {
-               n = ex_to<basic>(recombine_pair_to_ex(*it)).normal(sym_lst, repl_lst, level-1);
+               n = ex_to<basic>(recombine_pair_to_ex(*it)).normal(repl, level-1);
                num.push_back(n.op(0));
                den.push_back(n.op(1));
                it++;
        }
-       n = ex_to<numeric>(overall_coeff).normal(sym_lst, repl_lst, level-1);
+       n = ex_to<numeric>(overall_coeff).normal(repl, level-1);
        num.push_back(n.op(0));
        den.push_back(n.op(1));
 
@@ -1940,16 +1936,16 @@ ex mul::normal(lst &sym_lst, lst &repl_lst, int level) const
  *  distributes integer exponents to numerator and denominator, and replaces
  *  non-integer powers by temporary symbols.
  *  @see ex::normal */
-ex power::normal(lst &sym_lst, lst &repl_lst, int level) const
+ex power::normal(exmap & repl, int level) const
 {
        if (level == 1)
-               return (new lst(replace_with_symbol(*this, sym_lst, repl_lst), _ex1))->setflag(status_flags::dynallocated);
+               return (new lst(replace_with_symbol(*this, repl), _ex1))->setflag(status_flags::dynallocated);
        else if (level == -max_recursion_level)
                throw(std::runtime_error("max recursion level reached"));
 
        // Normalize basis and exponent (exponent gets reassembled)
-       ex n_basis = ex_to<basic>(basis).normal(sym_lst, repl_lst, level-1);
-       ex n_exponent = ex_to<basic>(exponent).normal(sym_lst, repl_lst, level-1);
+       ex n_basis = ex_to<basic>(basis).normal(repl, level-1);
+       ex n_exponent = ex_to<basic>(exponent).normal(repl, level-1);
        n_exponent = n_exponent.op(0) / n_exponent.op(1);
 
        if (n_exponent.info(info_flags::integer)) {
@@ -1970,32 +1966,32 @@ ex power::normal(lst &sym_lst, lst &repl_lst, int level) const
                if (n_exponent.info(info_flags::positive)) {
 
                        // (a/b)^x -> {sym((a/b)^x), 1}
-                       return (new lst(replace_with_symbol(power(n_basis.op(0) / n_basis.op(1), n_exponent), sym_lst, repl_lst), _ex1))->setflag(status_flags::dynallocated);
+                       return (new lst(replace_with_symbol(power(n_basis.op(0) / n_basis.op(1), n_exponent), repl), _ex1))->setflag(status_flags::dynallocated);
 
                } else if (n_exponent.info(info_flags::negative)) {
 
                        if (n_basis.op(1).is_equal(_ex1)) {
 
                                // a^-x -> {1, sym(a^x)}
-                               return (new lst(_ex1, replace_with_symbol(power(n_basis.op(0), -n_exponent), sym_lst, repl_lst)))->setflag(status_flags::dynallocated);
+                               return (new lst(_ex1, replace_with_symbol(power(n_basis.op(0), -n_exponent), repl)))->setflag(status_flags::dynallocated);
 
                        } else {
 
                                // (a/b)^-x -> {sym((b/a)^x), 1}
-                               return (new lst(replace_with_symbol(power(n_basis.op(1) / n_basis.op(0), -n_exponent), sym_lst, repl_lst), _ex1))->setflag(status_flags::dynallocated);
+                               return (new lst(replace_with_symbol(power(n_basis.op(1) / n_basis.op(0), -n_exponent), repl), _ex1))->setflag(status_flags::dynallocated);
                        }
                }
        }
 
        // (a/b)^x -> {sym((a/b)^x, 1}
-       return (new lst(replace_with_symbol(power(n_basis.op(0) / n_basis.op(1), n_exponent), sym_lst, repl_lst), _ex1))->setflag(status_flags::dynallocated);
+       return (new lst(replace_with_symbol(power(n_basis.op(0) / n_basis.op(1), n_exponent), repl), _ex1))->setflag(status_flags::dynallocated);
 }
 
 
 /** Implementation of ex::normal() for pseries. It normalizes each coefficient
  *  and replaces the series by a temporary symbol.
  *  @see ex::normal */
-ex pseries::normal(lst &sym_lst, lst &repl_lst, int level) const
+ex pseries::normal(exmap & repl, int level) const
 {
        epvector newseq;
        epvector::const_iterator i = seq.begin(), end = seq.end();
@@ -2006,7 +2002,7 @@ ex pseries::normal(lst &sym_lst, lst &repl_lst, int level) const
                ++i;
        }
        ex n = pseries(relational(var,point), newseq);
-       return (new lst(replace_with_symbol(n, sym_lst, repl_lst), _ex1))->setflag(status_flags::dynallocated);
+       return (new lst(replace_with_symbol(n, repl), _ex1))->setflag(status_flags::dynallocated);
 }
 
 
@@ -2024,14 +2020,14 @@ ex pseries::normal(lst &sym_lst, lst &repl_lst, int level) const
  *  @return normalized expression */
 ex ex::normal(int level) const
 {
-       lst sym_lst, repl_lst;
+       exmap repl;
 
-       ex e = bp->normal(sym_lst, repl_lst, level);
+       ex e = bp->normal(repl, level);
        GINAC_ASSERT(is_a<lst>(e));
 
        // Re-insert replaced symbols
-       if (sym_lst.nops() > 0)
-               e = e.subs(sym_lst, repl_lst);
+       if (!repl.empty())
+               e = e.subs(repl);
 
        // Convert {numerator, denominator} form back to fraction
        return e.op(0) / e.op(1);
@@ -2045,16 +2041,16 @@ ex ex::normal(int level) const
  *  @return numerator */
 ex ex::numer() const
 {
-       lst sym_lst, repl_lst;
+       exmap repl;
 
-       ex e = bp->normal(sym_lst, repl_lst, 0);
+       ex e = bp->normal(repl, 0);
        GINAC_ASSERT(is_a<lst>(e));
 
        // Re-insert replaced symbols
-       if (sym_lst.nops() > 0)
-               return e.op(0).subs(sym_lst, repl_lst);
-       else
+       if (repl.empty())
                return e.op(0);
+       else
+               return e.op(0).subs(repl);
 }
 
 /** Get denominator of an expression. If the expression is not of the normal
@@ -2065,16 +2061,16 @@ ex ex::numer() const
  *  @return denominator */
 ex ex::denom() const
 {
-       lst sym_lst, repl_lst;
+       exmap repl;
 
-       ex e = bp->normal(sym_lst, repl_lst, 0);
+       ex e = bp->normal(repl, 0);
        GINAC_ASSERT(is_a<lst>(e));
 
        // Re-insert replaced symbols
-       if (sym_lst.nops() > 0)
-               return e.op(1).subs(sym_lst, repl_lst);
-       else
+       if (repl.empty())
                return e.op(1);
+       else
+               return e.op(1).subs(repl);
 }
 
 /** Get numerator and denominator of an expression. If the expresison is not
@@ -2085,16 +2081,16 @@ ex ex::denom() const
  *  @return a list [numerator, denominator] */
 ex ex::numer_denom() const
 {
-       lst sym_lst, repl_lst;
+       exmap repl;
 
-       ex e = bp->normal(sym_lst, repl_lst, 0);
+       ex e = bp->normal(repl, 0);
        GINAC_ASSERT(is_a<lst>(e));
 
        // Re-insert replaced symbols
-       if (sym_lst.nops() > 0)
-               return e.subs(sym_lst, repl_lst);
-       else
+       if (repl.empty())
                return e;
+       else
+               return e.subs(repl);
 }
 
 
index 33498f3..a4f4447 100644 (file)
@@ -104,8 +104,8 @@ public:
        bool has(const ex &other) const;
        ex eval(int level = 0) const;
        ex evalf(int level = 0) const;
-       ex subs(const lst & ls, const lst & lr, unsigned options = 0) const { return subs_one_level(ls, lr, options); } // overwrites basic::subs() for performance reasons
-       ex normal(lst &sym_lst, lst &repl_lst, int level = 0) const;
+       ex subs(const exmap & m, unsigned options = 0) const { return subs_one_level(m, options); } // overwrites basic::subs() for performance reasons
+       ex normal(exmap & repl, int level = 0) const;
        ex to_rational(lst &repl_lst) const;
        ex to_polynomial(lst &repl_lst) const;
        numeric integer_content() const;
index 0d24a04..a30182e 100644 (file)
@@ -518,27 +518,26 @@ ex power::evalm() const
 // from mul.cpp
 extern bool tryfactsubs(const ex &, const ex &, int &, lst &);
 
-ex power::subs(const lst & ls, const lst & lr, unsigned options) const
+ex power::subs(const exmap & m, unsigned options) const
 {      
-       const ex &subsed_basis = basis.subs(ls, lr, options);
-       const ex &subsed_exponent = exponent.subs(ls, lr, options);
+       const ex &subsed_basis = basis.subs(m, options);
+       const ex &subsed_exponent = exponent.subs(m, options);
 
        if (!are_ex_trivially_equal(basis, subsed_basis)
         || !are_ex_trivially_equal(exponent, subsed_exponent)) 
-               return power(subsed_basis, subsed_exponent).subs_one_level(ls, lr, options);
+               return power(subsed_basis, subsed_exponent).subs_one_level(m, options);
 
        if (!(options & subs_options::subs_algebraic))
-               return subs_one_level(ls, lr, options);
+               return subs_one_level(m, options);
 
-       lst::const_iterator its, itr;
-       for (its = ls.begin(), itr = lr.begin(); its != ls.end(); ++its, ++itr) {
+       for (exmap::const_iterator it = m.begin(); it != m.end(); ++it) {
                int nummatches = std::numeric_limits<int>::max();
                lst repls;
-               if (tryfactsubs(*this, *its, nummatches, repls))
-                       return (ex_to<basic>((*this) * power(itr->subs(ex(repls), subs_options::subs_no_pattern) / its->subs(ex(repls), subs_options::subs_no_pattern), nummatches))).subs_one_level(ls, lr, options);
+               if (tryfactsubs(*this, it->first, nummatches, repls))
+                       return (ex_to<basic>((*this) * power(it->second.subs(ex(repls), subs_options::subs_no_pattern) / it->first.subs(ex(repls), subs_options::subs_no_pattern), nummatches))).subs_one_level(m, options);
        }
 
-       return subs_one_level(ls, lr, options);
+       return subs_one_level(m, options);
 }
 
 ex power::eval_ncmul(const exvector & v) const
index 4e1df9f..a1f5a23 100644 (file)
@@ -61,8 +61,8 @@ public:
        ex evalf(int level=0) const;
        ex evalm() const;
        ex series(const relational & s, int order, unsigned options = 0) const;
-       ex subs(const lst & ls, const lst & lr, unsigned options = 0) const;
-       ex normal(lst &sym_lst, lst &repl_lst, int level = 0) const;
+       ex subs(const exmap & m, unsigned options = 0) const;
+       ex normal(exmap & repl, int level = 0) const;
        ex to_rational(lst &repl_lst) const;
        ex to_polynomial(lst &repl_lst) const;
        exvector get_free_indices() const;
index e988705..22d7b77 100644 (file)
@@ -396,13 +396,13 @@ ex pseries::evalf(int level) const
        return (new pseries(relational(var,point), new_seq))->setflag(status_flags::dynallocated | status_flags::evaluated);
 }
 
-ex pseries::subs(const lst & ls, const lst & lr, unsigned options) const
+ex pseries::subs(const exmap & m, unsigned options) const
 {
        // If expansion variable is being substituted, convert the series to a
        // polynomial and do the substitution there because the result might
        // no longer be a power series
-       if (ls.has(var))
-               return convert_to_poly(true).subs(ls, lr, options);
+       if (m.find(var) != m.end())
+               return convert_to_poly(true).subs(m, options);
        
        // Otherwise construct a new series with substituted coefficients and
        // expansion point
@@ -410,10 +410,10 @@ ex pseries::subs(const lst & ls, const lst & lr, unsigned options) const
        newseq.reserve(seq.size());
        epvector::const_iterator it = seq.begin(), itend = seq.end();
        while (it != itend) {
-               newseq.push_back(expair(it->rest.subs(ls, lr, options), it->coeff));
+               newseq.push_back(expair(it->rest.subs(m, options), it->coeff));
                ++it;
        }
-       return (new pseries(relational(var,point.subs(ls, lr, options)), newseq))->setflag(status_flags::dynallocated);
+       return (new pseries(relational(var,point.subs(m, options)), newseq))->setflag(status_flags::dynallocated);
 }
 
 /** Implementation of ex::expand() for a power series.  It expands all the
index 17ed024..5dbdc66 100644 (file)
@@ -53,8 +53,8 @@ public:
        ex eval(int level=0) const;
        ex evalf(int level=0) const;
        ex series(const relational & r, int order, unsigned options = 0) const;
-       ex subs(const lst & ls, const lst & lr, unsigned options = 0) const;
-       ex normal(lst &sym_lst, lst &repl_lst, int level = 0) const;
+       ex subs(const exmap & m, unsigned options = 0) const;
+       ex normal(exmap & repl, int level = 0) const;
        ex expand(unsigned options = 0) const;
 protected:
        ex derivative(const symbol & s) const;
index c45dd1e..453b4c0 100644 (file)
@@ -178,15 +178,15 @@ ex relational::eval(int level) const
        return (new relational(lh.eval(level-1),rh.eval(level-1),o))->setflag(status_flags::dynallocated | status_flags::evaluated);
 }
 
-ex relational::subs(const lst & ls, const lst & lr, unsigned options) const
+ex relational::subs(const exmap & m, unsigned options) const
 {
-       const ex & subsed_lh = lh.subs(ls, lr, options);
-       const ex & subsed_rh = rh.subs(ls, lr, options);
+       const ex & subsed_lh = lh.subs(m, options);
+       const ex & subsed_rh = rh.subs(m, options);
 
        if (!are_ex_trivially_equal(lh, subsed_lh) || !are_ex_trivially_equal(rh, subsed_rh))
-               return relational(subsed_lh, subsed_rh, o).subs_one_level(ls, lr, options);
+               return relational(subsed_lh, subsed_rh, o).subs_one_level(m, options);
        else
-               return subs_one_level(ls, lr, options);
+               return subs_one_level(m, options);
 }
 
 ex relational::eval_ncmul(const exvector & v) const
index 591cc3b..cbe5146 100644 (file)
@@ -57,7 +57,7 @@ public:
        size_t nops() const;
        ex op(size_t i) const;
        ex map(map_function & f) const;
-       ex subs(const lst & ls, const lst & lr, unsigned options = 0) const;
+       ex subs(const exmap & m, unsigned options = 0) const;
        ex eval(int level=0) const;
 
 protected:
index 23ae93c..41b243a 100644 (file)
@@ -160,7 +160,7 @@ protected:
 public:
 
        // substitutions
-       ex subs(const lst & ls, const lst & lr, unsigned options = 0) const { return inherited::subs(ls, lr, options); }
+       ex subs(const exmap & m, unsigned options = 0) const { return inherited::subs(m, options); }
 
        // function mapping
        ex map(map_function & f) const { return inherited::map(f); }
@@ -181,7 +181,7 @@ public:
        ex series(const relational & r, int order, unsigned options = 0) const { return inherited::series(r, order, options); }
 
        // rational functions
-       ex normal(lst & sym_lst, lst & repl_lst, int level = 0) const { return inherited::normal(sym_lst, repl_lst, level); }
+       ex normal(exmap & repl, int level = 0) const { return inherited::normal(repl, level); }
        ex to_rational(lst & repl_lst) const { return inherited::to_rational(repl_lst); }
        ex to_polynomial(lst & repl_lst) const { return inherited::to_polynomial(repl_lst); }
 
index 8ff6313..8cd738b 100644 (file)
@@ -67,8 +67,8 @@ public:
        ex eval(int level = 0) const;
        ex evalf(int level = 0) const { return *this; } // overwrites basic::evalf() for performance reasons
        ex series(const relational & s, int order, unsigned options = 0) const;
-       ex subs(const lst & ls, const lst & lr, unsigned options = 0) const { return subs_one_level(ls, lr, options); } // overwrites basic::subs() for performance reasons
-       ex normal(lst &sym_lst, lst &repl_lst, int level = 0) const;
+       ex subs(const exmap & m, unsigned options = 0) const { return subs_one_level(m, options); } // overwrites basic::subs() for performance reasons
+       ex normal(exmap & repl, int level = 0) const;
        ex to_rational(lst &repl_lst) const;
        ex to_polynomial(lst &repl_lst) const;
 protected:
index d78eff6..f8e5635 100644 (file)
@@ -171,7 +171,10 @@ ex tensdelta::eval_indexed(const basic & i) const
        // dimension
        if (!i1.get_dim().is_equal(i2.get_dim())) {
                ex min_dim = i1.minimal_dim(i2);
-               return i.subs(lst(i1 == i1.replace_dim(min_dim), i2 == i2.replace_dim(min_dim)));
+               exmap m;
+               m[i1] = i1.replace_dim(min_dim);
+               m[i2] = i2.replace_dim(min_dim);
+               return i.subs(m);
        }
 
        // Trace of delta tensor is the (effective) dimension of the space
@@ -212,7 +215,10 @@ ex tensmetric::eval_indexed(const basic & i) const
        // dimension
        if (!i1.get_dim().is_equal(i2.get_dim())) {
                ex min_dim = i1.minimal_dim(i2);
-               return i.subs(lst(i1 == i1.replace_dim(min_dim), i2 == i2.replace_dim(min_dim)));
+               exmap m;
+               m[i1] = i1.replace_dim(min_dim);
+               m[i2] = i2.replace_dim(min_dim);
+               return i.subs(m);
        }
 
        // A metric tensor with one covariant and one contravariant index gets