]> www.ginac.de Git - ginac.git/blobdiff - ginac/normal.cpp
Convert some more ctors from copying to moving STL containers.
[ginac.git] / ginac / normal.cpp
index e97a7ce4a2c19c7f530881e5c1257e7562f20604..e30afda6b941b90b4f1ee97aa59595d2e4bb9c0a 100644 (file)
@@ -6,7 +6,7 @@
  *  computation, square-free factorization and rational function normalization. */
 
 /*
- *  GiNaC Copyright (C) 1999-2008 Johannes Gutenberg University Mainz, Germany
+ *  GiNaC Copyright (C) 1999-2015 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
@@ -23,9 +23,6 @@
  *  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
  */
 
-#include <algorithm>
-#include <map>
-
 #include "normal.h"
 #include "basic.h"
 #include "ex.h"
 #include "pseries.h"
 #include "symbol.h"
 #include "utils.h"
+#include "polynomial/chinrem_gcd.h"
+
+#include <algorithm>
+#include <map>
 
 namespace GiNaC {
 
@@ -119,6 +120,11 @@ static bool get_first_symbol(const ex &e, ex &x)
  *
  *  @see get_symbol_stats */
 struct sym_desc {
+       /** Initialize symbol, leave other variables uninitialized */
+       sym_desc(const ex& s)
+         : sym(s), deg_a(0), deg_b(0), ldeg_a(0), ldeg_b(0), max_deg(0), max_lcnops(0)
+       { }
+
        /** Reference to symbol */
        ex sym;
 
@@ -156,15 +162,11 @@ typedef std::vector<sym_desc> sym_desc_vec;
 // Add symbol the sym_desc_vec (used internally by get_symbol_stats())
 static void add_symbol(const ex &s, sym_desc_vec &v)
 {
-       sym_desc_vec::const_iterator it = v.begin(), itend = v.end();
-       while (it != itend) {
-               if (it->sym.is_equal(s))  // If it's already in there, don't add it a second time
+       for (auto & it : v)
+               if (it.sym.is_equal(s))  // If it's already in there, don't add it a second time
                        return;
-               ++it;
-       }
-       sym_desc d;
-       d.sym = s;
-       v.push_back(d);
+
+       v.push_back(sym_desc(s));
 }
 
 // Collect all symbols of an expression (used internally by get_symbol_stats())
@@ -196,17 +198,15 @@ static void get_symbol_stats(const ex &a, const ex &b, sym_desc_vec &v)
 {
        collect_symbols(a.eval(), v);   // eval() to expand assigned symbols
        collect_symbols(b.eval(), v);
-       sym_desc_vec::iterator it = v.begin(), itend = v.end();
-       while (it != itend) {
-               int deg_a = a.degree(it->sym);
-               int deg_b = b.degree(it->sym);
-               it->deg_a = deg_a;
-               it->deg_b = deg_b;
-               it->max_deg = std::max(deg_a, deg_b);
-               it->max_lcnops = std::max(a.lcoeff(it->sym).nops(), b.lcoeff(it->sym).nops());
-               it->ldeg_a = a.ldegree(it->sym);
-               it->ldeg_b = b.ldegree(it->sym);
-               ++it;
+       for (auto & it : v) {
+               int deg_a = a.degree(it.sym);
+               int deg_b = b.degree(it.sym);
+               it.deg_a = deg_a;
+               it.deg_b = deg_b;
+               it.max_deg = std::max(deg_a, deg_b);
+               it.max_lcnops = std::max(a.lcoeff(it.sym).nops(), b.lcoeff(it.sym).nops());
+               it.ldeg_a = a.ldegree(it.sym);
+               it.ldeg_b = b.ldegree(it.sym);
        }
        std::sort(v.begin(), v.end());
 
@@ -320,15 +320,12 @@ numeric numeric::integer_content() const
 
 numeric add::integer_content() const
 {
-       epvector::const_iterator it = seq.begin();
-       epvector::const_iterator itend = seq.end();
        numeric c = *_num0_p, l = *_num1_p;
-       while (it != itend) {
-               GINAC_ASSERT(!is_exactly_a<numeric>(it->rest));
-               GINAC_ASSERT(is_exactly_a<numeric>(it->coeff));
-               c = gcd(ex_to<numeric>(it->coeff).numer(), c);
-               l = lcm(ex_to<numeric>(it->coeff).denom(), l);
-               it++;
+       for (auto & it : seq) {
+               GINAC_ASSERT(!is_exactly_a<numeric>(it.rest));
+               GINAC_ASSERT(is_exactly_a<numeric>(it.coeff));
+               c = gcd(ex_to<numeric>(it.coeff).numer(), c);
+               l = lcm(ex_to<numeric>(it.coeff).denom(), l);
        }
        GINAC_ASSERT(is_exactly_a<numeric>(overall_coeff));
        c = gcd(ex_to<numeric>(overall_coeff).numer(), c);
@@ -339,11 +336,8 @@ numeric add::integer_content() const
 numeric mul::integer_content() const
 {
 #ifdef DO_GINAC_ASSERT
-       epvector::const_iterator it = seq.begin();
-       epvector::const_iterator itend = seq.end();
-       while (it != itend) {
-               GINAC_ASSERT(!is_exactly_a<numeric>(recombine_pair_to_ex(*it)));
-               ++it;
+       for (auto & it : seq) {
+               GINAC_ASSERT(!is_exactly_a<numeric>(recombine_pair_to_ex(it)));
        }
 #endif // def DO_GINAC_ASSERT
        GINAC_ASSERT(is_exactly_a<numeric>(overall_coeff));
@@ -672,12 +666,13 @@ bool divide(const ex &a, const ex &b, ex &q, bool check_args)
                        q = rem_i*power(ab, a_exp - 1);
                        return true;
                }
-               for (int i=2; i < a_exp; i++) {
-                       if (divide(power(ab, i), b, rem_i, false)) {
-                               q = rem_i*power(ab, a_exp - i);
-                               return true;
-                       }
-               } // ... so we *really* need to expand expression.
+// code below is commented-out because it leads to a significant slowdown
+//             for (int i=2; i < a_exp; i++) {
+//                     if (divide(power(ab, i), b, rem_i, false)) {
+//                             q = rem_i*power(ab, a_exp - i);
+//                             return true;
+//                     }
+//             } // ... so we *really* need to expand expression.
        }
        
        // Polynomial long division (recursive)
@@ -794,10 +789,10 @@ static bool divide_in_z(const ex &a, const ex &b, ex &q, sym_desc_vec::const_ite
 
        if (is_exactly_a<mul>(b)) {
                ex qbar = a;
-               for (const_iterator itrb = b.begin(); itrb != b.end(); ++itrb) {
+               for (const auto & it : b) {
                        sym_desc_vec sym_stats;
-                       get_symbol_stats(a, *itrb, sym_stats);
-                       if (!divide_in_z(qbar, *itrb, q, sym_stats.begin()))
+                       get_symbol_stats(a, it, sym_stats);
+                       if (!divide_in_z(qbar, it, q, sym_stats.begin()))
                                return false;
 
                        qbar = q;
@@ -959,7 +954,7 @@ ex ex::content(const ex &x) const
                return lcoeff * c / lcoeff.unit(x);
        ex cont = _ex0;
        for (int i=ldeg; i<=deg; i++)
-               cont = gcd(r.coeff(x, i), cont, NULL, NULL, false);
+               cont = gcd(r.coeff(x, i), cont, nullptr, nullptr, false);
        return cont * c;
 }
 
@@ -1099,7 +1094,7 @@ static ex sr_gcd(const ex &a, const ex &b, sym_desc_vec::const_iterator var)
        // Remove content from c and d, to be attached to GCD later
        ex cont_c = c.content(x);
        ex cont_d = d.content(x);
-       ex gamma = gcd(cont_c, cont_d, NULL, NULL, false);
+       ex gamma = gcd(cont_c, cont_d, nullptr, nullptr, false);
        if (ddeg == 0)
                return gamma;
        c = c.primpart(x, cont_c);
@@ -1163,17 +1158,14 @@ numeric numeric::max_coefficient() const
 
 numeric add::max_coefficient() const
 {
-       epvector::const_iterator it = seq.begin();
-       epvector::const_iterator itend = seq.end();
        GINAC_ASSERT(is_exactly_a<numeric>(overall_coeff));
        numeric cur_max = abs(ex_to<numeric>(overall_coeff));
-       while (it != itend) {
+       for (auto & it : seq) {
                numeric a;
-               GINAC_ASSERT(!is_exactly_a<numeric>(it->rest));
-               a = abs(ex_to<numeric>(it->coeff));
+               GINAC_ASSERT(!is_exactly_a<numeric>(it.rest));
+               a = abs(ex_to<numeric>(it.coeff));
                if (a > cur_max)
                        cur_max = a;
-               it++;
        }
        return cur_max;
 }
@@ -1181,11 +1173,8 @@ numeric add::max_coefficient() const
 numeric mul::max_coefficient() const
 {
 #ifdef DO_GINAC_ASSERT
-       epvector::const_iterator it = seq.begin();
-       epvector::const_iterator itend = seq.end();
-       while (it != itend) {
-               GINAC_ASSERT(!is_exactly_a<numeric>(recombine_pair_to_ex(*it)));
-               it++;
+       for (auto & it : seq) {
+               GINAC_ASSERT(!is_exactly_a<numeric>(recombine_pair_to_ex(it)));
        }
 #endif // def DO_GINAC_ASSERT
        GINAC_ASSERT(is_exactly_a<numeric>(overall_coeff));
@@ -1213,28 +1202,22 @@ ex add::smod(const numeric &xi) const
 {
        epvector newseq;
        newseq.reserve(seq.size()+1);
-       epvector::const_iterator it = seq.begin();
-       epvector::const_iterator itend = seq.end();
-       while (it != itend) {
-               GINAC_ASSERT(!is_exactly_a<numeric>(it->rest));
-               numeric coeff = GiNaC::smod(ex_to<numeric>(it->coeff), xi);
+       for (auto & it : seq) {
+               GINAC_ASSERT(!is_exactly_a<numeric>(it.rest));
+               numeric coeff = GiNaC::smod(ex_to<numeric>(it.coeff), xi);
                if (!coeff.is_zero())
-                       newseq.push_back(expair(it->rest, coeff));
-               it++;
+                       newseq.push_back(expair(it.rest, coeff));
        }
        GINAC_ASSERT(is_exactly_a<numeric>(overall_coeff));
        numeric coeff = GiNaC::smod(ex_to<numeric>(overall_coeff), xi);
-       return (new add(newseq,coeff))->setflag(status_flags::dynallocated);
+       return (new add(std::move(newseq), coeff))->setflag(status_flags::dynallocated);
 }
 
 ex mul::smod(const numeric &xi) const
 {
 #ifdef DO_GINAC_ASSERT
-       epvector::const_iterator it = seq.begin();
-       epvector::const_iterator itend = seq.end();
-       while (it != itend) {
-               GINAC_ASSERT(!is_exactly_a<numeric>(recombine_pair_to_ex(*it)));
-               it++;
+       for (auto & it : seq) {
+               GINAC_ASSERT(!is_exactly_a<numeric>(recombine_pair_to_ex(it)));
        }
 #endif // def DO_GINAC_ASSERT
        mul * mulcopyp = new mul(*this);
@@ -1270,9 +1253,9 @@ class gcdheu_failed {};
  *
  *  @param a  first integer multivariate polynomial (expanded)
  *  @param b  second integer multivariate polynomial (expanded)
- *  @param ca  cofactor of polynomial a (returned), NULL to suppress
+ *  @param ca  cofactor of polynomial a (returned), nullptr to suppress
  *             calculation of cofactor
- *  @param cb  cofactor of polynomial b (returned), NULL to suppress
+ *  @param cb  cofactor of polynomial b (returned), nullptr to suppress
  *             calculation of cofactor
  *  @param var iterator to first element of vector of sym_desc structs
  *  @param res the GCD (returned)
@@ -1363,9 +1346,9 @@ static bool heur_gcd_z(ex& res, const ex &a, const ex &b, ex *ca, ex *cb,
  *
  *  @param a  first rational multivariate polynomial (expanded)
  *  @param b  second rational multivariate polynomial (expanded)
- *  @param ca  cofactor of polynomial a (returned), NULL to suppress
+ *  @param ca  cofactor of polynomial a (returned), nullptr to suppress
  *             calculation of cofactor
- *  @param cb  cofactor of polynomial b (returned), NULL to suppress
+ *  @param cb  cofactor of polynomial b (returned), nullptr to suppress
  *             calculation of cofactor
  *  @param var iterator to first element of vector of sym_desc structs
  *  @param res the GCD (returned)
@@ -1429,8 +1412,8 @@ static ex gcd_pf_mul(const ex& a, const ex& b, ex* ca, ex* cb);
  *
  *  @param a  first multivariate polynomial
  *  @param b  second multivariate polynomial
- *  @param ca pointer to expression that will receive the cofactor of a, or NULL
- *  @param cb pointer to expression that will receive the cofactor of b, or NULL
+ *  @param ca pointer to expression that will receive the cofactor of a, or nullptr
+ *  @param cb pointer to expression that will receive the cofactor of b, or nullptr
  *  @param check_args  check whether a and b are polynomials with rational
  *         coefficients (defaults to "true")
  *  @return the GCD as a new expression */
@@ -1533,7 +1516,7 @@ ex gcd(const ex &a, const ex &b, ex *ca, ex *cb, bool check_args, unsigned optio
                if (ca)
                        *ca = ex_to<numeric>(aex)/g;
                if (cb)
-                       *cb = bex/g;
+                       *cb = bex/g;
                return g;
        }
 
@@ -1622,8 +1605,15 @@ ex gcd(const ex &a, const ex &b, ex *ca, ex *cb, bool check_args, unsigned optio
                }
 #endif
        }
+       if (options & gcd_options::use_sr_gcd) {
+               g = sr_gcd(aex, bex, var);
+       } else {
+               exvector vars;
+               for (std::size_t n = sym_stats.size(); n-- != 0; )
+                       vars.push_back(sym_stats[n].sym);
+               g = chinrem_gcd(aex, bex, vars);
+       }
 
-       g = sr_gcd(aex, bex, var);
        if (g.is_equal(_ex1)) {
                // Keep cofactors factored if possible
                if (ca)
@@ -1639,111 +1629,93 @@ ex gcd(const ex &a, const ex &b, ex *ca, ex *cb, bool check_args, unsigned optio
        return g;
 }
 
-static ex gcd_pf_pow(const ex& a, const ex& b, ex* ca, ex* cb)
+// gcd helper to handle partially factored polynomials (to avoid expanding
+// large expressions). Both arguments should be powers.
+static ex gcd_pf_pow_pow(const ex& a, const ex& b, ex* ca, ex* cb)
 {
-       if (is_exactly_a<power>(a)) {
-               ex p = a.op(0);
-               const ex& exp_a = a.op(1);
-               if (is_exactly_a<power>(b)) {
-                       ex pb = b.op(0);
-                       const ex& exp_b = b.op(1);
-                       if (p.is_equal(pb)) {
-                               // a = p^n, b = p^m, gcd = p^min(n, m)
-                               if (exp_a < exp_b) {
-                                       if (ca)
-                                               *ca = _ex1;
-                                       if (cb)
-                                               *cb = power(p, exp_b - exp_a);
-                                       return power(p, exp_a);
-                               } else {
-                                       if (ca)
-                                               *ca = power(p, exp_a - exp_b);
-                                       if (cb)
-                                               *cb = _ex1;
-                                       return power(p, exp_b);
-                               }
-                       } else {
-                               ex p_co, pb_co;
-                               ex p_gcd = gcd(p, pb, &p_co, &pb_co, false);
-                               if (p_gcd.is_equal(_ex1)) {
-                                       // a(x) = p(x)^n, b(x) = p_b(x)^m, gcd (p, p_b) = 1 ==>
-                                       // gcd(a,b) = 1
-                                       if (ca)
-                                               *ca = a;
-                                       if (cb)
-                                               *cb = b;
-                                       return _ex1;
-                                       // XXX: do I need to check for p_gcd = -1?
-                               } else {
-                                       // there are common factors:
-                                       // a(x) = g(x)^n A(x)^n, b(x) = g(x)^m B(x)^m ==>
-                                       // gcd(a, b) = g(x)^n gcd(A(x)^n, g(x)^(n-m) B(x)^m
-                                       if (exp_a < exp_b) {
-                                               return power(p_gcd, exp_a)*
-                                                       gcd(power(p_co, exp_a), power(p_gcd, exp_b-exp_a)*power(pb_co, exp_b), ca, cb, false);
-                                       } else {
-                                               return power(p_gcd, exp_b)*
-                                                       gcd(power(p_gcd, exp_a - exp_b)*power(p_co, exp_a), power(pb_co, exp_b), ca, cb, false);
-                                       }
-                               } // p_gcd.is_equal(_ex1)
-                       } // p.is_equal(pb)
-
-               } else {
-                       if (p.is_equal(b)) {
-                               // a = p^n, b = p, gcd = p
-                               if (ca)
-                                       *ca = power(p, a.op(1) - 1);
-                               if (cb)
-                                       *cb = _ex1;
-                               return p;
-                       } 
+       ex p = a.op(0);
+       const ex& exp_a = a.op(1);
+       ex pb = b.op(0);
+       const ex& exp_b = b.op(1);
 
-                       ex p_co, bpart_co;
-                       ex p_gcd = gcd(p, b, &p_co, &bpart_co, false);
-
-                       if (p_gcd.is_equal(_ex1)) {
-                               // a(x) = p(x)^n, gcd(p, b) = 1 ==> gcd(a, b) = 1
-                               if (ca)
-                                       *ca = a;
-                               if (cb)
-                                       *cb = b;
-                               return _ex1;
-                       } else {
-                               // a(x) = g(x)^n A(x)^n, b(x) = g(x) B(x) ==> gcd(a, b) = g(x) gcd(g(x)^(n-1) A(x)^n, B(x))
-                               return p_gcd*gcd(power(p_gcd, exp_a-1)*power(p_co, exp_a), bpart_co, ca, cb, false);
-                       }
-               } // is_exactly_a<power>(b)
-
-       } else if (is_exactly_a<power>(b)) {
-               ex p = b.op(0);
-               if (p.is_equal(a)) {
-                       // a = p, b = p^n, gcd = p
+       // a = p^n, b = p^m, gcd = p^min(n, m)
+       if (p.is_equal(pb)) {
+               if (exp_a < exp_b) {
                        if (ca)
                                *ca = _ex1;
                        if (cb)
-                               *cb = power(p, b.op(1) - 1);
-                       return p;
+                               *cb = power(p, exp_b - exp_a);
+                       return power(p, exp_a);
+               } else {
+                       if (ca)
+                               *ca = power(p, exp_a - exp_b);
+                       if (cb)
+                               *cb = _ex1;
+                       return power(p, exp_b);
                }
+       }
 
-               ex p_co, apart_co;
-               const ex& exp_b(b.op(1));
-               ex p_gcd = gcd(a, p, &apart_co, &p_co, false);
-               if (p_gcd.is_equal(_ex1)) {
-                       // b=p(x)^n, gcd(a, p) = 1 ==> gcd(a, b) == 1
+       ex p_co, pb_co;
+       ex p_gcd = gcd(p, pb, &p_co, &pb_co, false);
+       // a(x) = p(x)^n, b(x) = p_b(x)^m, gcd (p, p_b) = 1 ==> gcd(a,b) = 1
+       if (p_gcd.is_equal(_ex1)) {
                        if (ca)
                                *ca = a;
                        if (cb)
                                *cb = b;
                        return _ex1;
-               } else {
-                       // there are common factors:
-                       // a(x) = g(x) A(x), b(x) = g(x)^n B(x)^n ==> gcd = g(x) gcd(g(x)^(n-1) A(x)^n, B(x))
+                       // XXX: do I need to check for p_gcd = -1?
+       }
 
-                       return p_gcd*gcd(apart_co, power(p_gcd, exp_b-1)*power(p_co, exp_b), ca, cb, false);
-               } // p_gcd.is_equal(_ex1)
+       // there are common factors:
+       // a(x) = g(x)^n A(x)^n, b(x) = g(x)^m B(x)^m ==>
+       // gcd(a, b) = g(x)^n gcd(A(x)^n, g(x)^(n-m) B(x)^m
+       if (exp_a < exp_b) {
+               ex pg =  gcd(power(p_co, exp_a), power(p_gcd, exp_b-exp_a)*power(pb_co, exp_b), ca, cb, false);
+               return power(p_gcd, exp_a)*pg;
+       } else {
+               ex pg = gcd(power(p_gcd, exp_a - exp_b)*power(p_co, exp_a), power(pb_co, exp_b), ca, cb, false);
+               return power(p_gcd, exp_b)*pg;
        }
 }
 
+static ex gcd_pf_pow(const ex& a, const ex& b, ex* ca, ex* cb)
+{
+       if (is_exactly_a<power>(a) && is_exactly_a<power>(b))
+               return gcd_pf_pow_pow(a, b, ca, cb);
+
+       if (is_exactly_a<power>(b) && (! is_exactly_a<power>(a)))
+               return gcd_pf_pow(b, a, cb, ca);
+
+       GINAC_ASSERT(is_exactly_a<power>(a));
+
+       ex p = a.op(0);
+       const ex& exp_a = a.op(1);
+       if (p.is_equal(b)) {
+               // a = p^n, b = p, gcd = p
+               if (ca)
+                       *ca = power(p, a.op(1) - 1);
+               if (cb)
+                       *cb = _ex1;
+               return p;
+       } 
+
+       ex p_co, bpart_co;
+       ex p_gcd = gcd(p, b, &p_co, &bpart_co, false);
+
+       // a(x) = p(x)^n, gcd(p, b) = 1 ==> gcd(a, b) = 1
+       if (p_gcd.is_equal(_ex1)) {
+               if (ca)
+                       *ca = a;
+               if (cb)
+                       *cb = b;
+               return _ex1;
+       }
+       // a(x) = g(x)^n A(x)^n, b(x) = g(x) B(x) ==> gcd(a, b) = g(x) gcd(g(x)^(n-1) A(x)^n, B(x))
+       ex rg = gcd(power(p_gcd, exp_a-1)*power(p_co, exp_a), bpart_co, ca, cb, false);
+       return p_gcd*rg;
+}
+
 static ex gcd_pf_mul(const ex& a, const ex& b, ex* ca, ex* cb)
 {
        if (is_exactly_a<mul>(a) && is_exactly_a<mul>(b)
@@ -1872,11 +1844,8 @@ ex sqrfree(const ex &a, const lst &l)
        if (l.nops()==0) {
                sym_desc_vec sdv;
                get_symbol_stats(a, _ex0, sdv);
-               sym_desc_vec::const_iterator it = sdv.begin(), itend = sdv.end();
-               while (it != itend) {
-                       args.append(it->sym);
-                       ++it;
-               }
+               for (auto & it : sdv)
+                       args.append(it.sym);
        } else {
                args = l;
        }
@@ -1899,18 +1868,15 @@ ex sqrfree(const ex &a, const lst &l)
 
        // recurse down the factors in remaining variables
        if (newargs.nops()>0) {
-               exvector::iterator i = factors.begin();
-               while (i != factors.end()) {
-                       *i = sqrfree(*i, newargs);
-                       ++i;
-               }
+               for (auto & it : factors)
+                       it = sqrfree(it, newargs);
        }
 
        // Done with recursion, now construct the final result
        ex result = _ex1;
-       exvector::const_iterator it = factors.begin(), itend = factors.end();
-       for (int p = 1; it!=itend; ++it, ++p)
-               result *= power(*it, p);
+       int p = 1;
+       for (auto & it : factors)
+               result *= power(it, p++);
 
        // Yun's algorithm does not account for constant factors.  (For univariate
        // polynomials it works only in the monic case.)  We can correct this by
@@ -1921,7 +1887,7 @@ ex sqrfree(const ex &a, const lst &l)
        else
                result *= quo(tmp, result, x);
 
-       // Put in the reational overall factor again and return
+       // Put in the rational overall factor again and return
        return result * lcm.inverse();
 }
 
@@ -2015,16 +1981,18 @@ ex sqrfree_parfrac(const ex & a, const symbol & x)
  *  @see ex::normal */
 static ex replace_with_symbol(const ex & e, exmap & repl, exmap & rev_lookup)
 {
+       // Since the repl contains replaced expressions we should search for them
+       ex e_replaced = e.subs(repl, subs_options::no_pattern);
+
        // Expression already replaced? Then return the assigned symbol
-       exmap::const_iterator it = rev_lookup.find(e);
+       auto it = rev_lookup.find(e_replaced);
        if (it != rev_lookup.end())
                return it->second;
-       
+
        // Otherwise create new symbol and add to list, taking care that the
        // replacement expression doesn't itself contain symbols from repl,
        // because subs() is not recursive
        ex es = (new symbol)->setflag(status_flags::dynallocated);
-       ex e_replaced = e.subs(repl, subs_options::no_pattern);
        repl.insert(std::make_pair(es, e_replaced));
        rev_lookup.insert(std::make_pair(e_replaced, es));
        return es;
@@ -2037,16 +2005,18 @@ static ex replace_with_symbol(const ex & e, exmap & repl, exmap & rev_lookup)
  *  @see basic::to_polynomial */
 static ex replace_with_symbol(const ex & e, exmap & repl)
 {
+       // Since the repl contains replaced expressions we should search for them
+       ex e_replaced = e.subs(repl, subs_options::no_pattern);
+
        // Expression already replaced? 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;
-       
+       for (auto & it : repl)
+               if (it.second.is_equal(e_replaced))
+                       return it.first;
+
        // Otherwise create new symbol and add to list, taking care that the
        // replacement expression doesn't itself contain symbols from repl,
        // because subs() is not recursive
        ex es = (new symbol)->setflag(status_flags::dynallocated);
-       ex e_replaced = e.subs(repl, subs_options::no_pattern);
        repl.insert(std::make_pair(es, e_replaced));
        return es;
 }
@@ -2186,12 +2156,10 @@ ex add::normal(exmap & repl, exmap & rev_lookup, int level) const
        exvector nums, dens;
        nums.reserve(seq.size()+1);
        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(repl, rev_lookup, level-1);
+       for (auto & it : seq) {
+               ex n = ex_to<basic>(recombine_pair_to_ex(it)).normal(repl, rev_lookup, level-1);
                nums.push_back(n.op(0));
                dens.push_back(n.op(1));
-               it++;
        }
        ex n = ex_to<numeric>(overall_coeff).normal(repl, rev_lookup, level-1);
        nums.push_back(n.op(0));
@@ -2203,8 +2171,8 @@ ex add::normal(exmap & repl, exmap & rev_lookup, int level) const
 //std::clog << "add::normal uses " << nums.size() << " summands:\n";
 
        // Add fractions sequentially
-       exvector::const_iterator num_it = nums.begin(), num_itend = nums.end();
-       exvector::const_iterator den_it = dens.begin(), den_itend = dens.end();
+       auto num_it = nums.begin(), num_itend = nums.end();
+       auto den_it = dens.begin(), den_itend = dens.end();
 //std::clog << " num = " << *num_it << ", den = " << *den_it << std::endl;
        ex num = *num_it++, den = *den_it++;
        while (num_it != num_itend) {
@@ -2217,7 +2185,7 @@ ex add::normal(exmap & repl, exmap & rev_lookup, int level) const
                        num_it++; den_it++;
                }
 
-               // Additiion of two fractions, taking advantage of the fact that
+               // Addition of two fractions, taking advantage of the fact that
                // the heuristic GCD algorithm computes the cofactors at no extra cost
                ex co_den1, co_den2;
                ex g = gcd(den, next_den, &co_den1, &co_den2, false);
@@ -2245,12 +2213,10 @@ ex mul::normal(exmap & repl, exmap & rev_lookup, int level) const
        exvector num; num.reserve(seq.size());
        exvector den; den.reserve(seq.size());
        ex n;
-       epvector::const_iterator it = seq.begin(), itend = seq.end();
-       while (it != itend) {
-               n = ex_to<basic>(recombine_pair_to_ex(*it)).normal(repl, rev_lookup, level-1);
+       for (auto & it : seq) {
+               n = ex_to<basic>(recombine_pair_to_ex(it)).normal(repl, rev_lookup, level-1);
                num.push_back(n.op(0));
                den.push_back(n.op(1));
-               it++;
        }
        n = ex_to<numeric>(overall_coeff).normal(repl, rev_lookup, level-1);
        num.push_back(n.op(0));
@@ -2324,14 +2290,12 @@ ex power::normal(exmap & repl, exmap & rev_lookup, int level) const
 ex pseries::normal(exmap & repl, exmap & rev_lookup, int level) const
 {
        epvector newseq;
-       epvector::const_iterator i = seq.begin(), end = seq.end();
-       while (i != end) {
-               ex restexp = i->rest.normal();
+       for (auto & it : seq) {
+               ex restexp = it.rest.normal();
                if (!restexp.is_zero())
-                       newseq.push_back(expair(restexp, i->coeff));
-               ++i;
+                       newseq.push_back(expair(restexp, it.coeff));
        }
-       ex n = pseries(relational(var,point), newseq);
+       ex n = pseries(relational(var,point), std::move(newseq));
        return (new lst(replace_with_symbol(n, repl, rev_lookup), _ex1))->setflag(status_flags::dynallocated);
 }
 
@@ -2403,7 +2367,7 @@ ex ex::denom() const
                return e.op(1).subs(repl, subs_options::no_pattern);
 }
 
-/** Get numerator and denominator of an expression. If the expresison is not
+/** Get numerator and denominator of an expression. If the expression is not
  *  of the normal form "numerator/denominator", it is first converted to this
  *  form and then a list [numerator, denominator] is returned.
  *
@@ -2447,15 +2411,15 @@ ex ex::to_rational(lst & repl_lst) const
 {
        // Convert lst to exmap
        exmap m;
-       for (lst::const_iterator it = repl_lst.begin(); it != repl_lst.end(); ++it)
-               m.insert(std::make_pair(it->op(0), it->op(1)));
+       for (auto & it : repl_lst)
+               m.insert(std::make_pair(it.op(0), it.op(1)));
 
        ex ret = bp->to_rational(m);
 
        // Convert exmap back to lst
        repl_lst.remove_all();
-       for (exmap::const_iterator it = m.begin(); it != m.end(); ++it)
-               repl_lst.append(it->first == it->second);
+       for (auto & it : m)
+               repl_lst.append(it.first == it.second);
 
        return ret;
 }
@@ -2470,15 +2434,15 @@ ex ex::to_polynomial(lst & repl_lst) const
 {
        // Convert lst to exmap
        exmap m;
-       for (lst::const_iterator it = repl_lst.begin(); it != repl_lst.end(); ++it)
-               m.insert(std::make_pair(it->op(0), it->op(1)));
+       for (auto & it : repl_lst)
+               m.insert(std::make_pair(it.op(0), it.op(1)));
 
        ex ret = bp->to_polynomial(m);
 
        // Convert exmap back to lst
        repl_lst.remove_all();
-       for (exmap::const_iterator it = m.begin(); it != m.end(); ++it)
-               repl_lst.append(it->first == it->second);
+       for (auto & it : m)
+               repl_lst.append(it.first == it.second);
 
        return ret;
 }
@@ -2585,17 +2549,15 @@ ex expairseq::to_rational(exmap & repl) const
 {
        epvector s;
        s.reserve(seq.size());
-       epvector::const_iterator i = seq.begin(), end = seq.end();
-       while (i != end) {
-               s.push_back(split_ex_to_pair(recombine_pair_to_ex(*i).to_rational(repl)));
-               ++i;
-       }
+       for (auto & it : seq)
+               s.push_back(split_ex_to_pair(recombine_pair_to_ex(it).to_rational(repl)));
+
        ex oc = overall_coeff.to_rational(repl);
        if (oc.info(info_flags::numeric))
-               return thisexpairseq(s, overall_coeff);
+               return thisexpairseq(std::move(s), overall_coeff);
        else
                s.push_back(combine_ex_with_coeff_to_pair(oc, _ex1));
-       return thisexpairseq(s, default_overall_coeff());
+       return thisexpairseq(std::move(s), default_overall_coeff());
 }
 
 /** Implementation of ex::to_polynomial() for expairseqs. */
@@ -2603,17 +2565,15 @@ ex expairseq::to_polynomial(exmap & repl) const
 {
        epvector s;
        s.reserve(seq.size());
-       epvector::const_iterator i = seq.begin(), end = seq.end();
-       while (i != end) {
-               s.push_back(split_ex_to_pair(recombine_pair_to_ex(*i).to_polynomial(repl)));
-               ++i;
-       }
+       for (auto & it : seq)
+               s.push_back(split_ex_to_pair(recombine_pair_to_ex(it).to_polynomial(repl)));
+
        ex oc = overall_coeff.to_polynomial(repl);
        if (oc.info(info_flags::numeric))
-               return thisexpairseq(s, overall_coeff);
+               return thisexpairseq(std::move(s), overall_coeff);
        else
                s.push_back(combine_ex_with_coeff_to_pair(oc, _ex1));
-       return thisexpairseq(s, default_overall_coeff());
+       return thisexpairseq(std::move(s), default_overall_coeff());
 }