]> www.ginac.de Git - ginac.git/blobdiff - ginac/indexed.cpp
Improve method of setting status_flags::dynallocated.
[ginac.git] / ginac / indexed.cpp
index d772b3f4743bd85ff8c127110f888c34801c0c6c..ca4beff53f63d2b57761c4380ae68a9b89829b1c 100644 (file)
@@ -3,7 +3,7 @@
  *  Implementation of GiNaC's indexed expressions. */
 
 /*
- *  GiNaC Copyright (C) 1999-2006 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
  *  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
  */
 
-#include <iostream>
-#include <sstream>
-#include <stdexcept>
-
 #include "indexed.h"
 #include "idx.h"
 #include "add.h"
 #include "matrix.h"
 #include "inifcns.h"
 
+#include <iostream>
+#include <limits>
+#include <sstream>
+#include <stdexcept>
+
 namespace GiNaC {
 
 GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(indexed, exprseq,
@@ -54,96 +55,83 @@ GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(indexed, exprseq,
 
 indexed::indexed() : symtree(not_symmetric())
 {
-       tinfo_key = &indexed::tinfo_static;
 }
 
 //////////
 // other constructors
 //////////
 
-indexed::indexed(const ex & b) : inherited(b), symtree(not_symmetric())
+indexed::indexed(const ex & b) : inherited{b}, symtree(not_symmetric())
 {
-       tinfo_key = &indexed::tinfo_static;
        validate();
 }
 
-indexed::indexed(const ex & b, const ex & i1) : inherited(b, i1), symtree(not_symmetric())
+indexed::indexed(const ex & b, const ex & i1) : inherited{b, i1}, symtree(not_symmetric())
 {
-       tinfo_key = &indexed::tinfo_static;
        validate();
 }
 
-indexed::indexed(const ex & b, const ex & i1, const ex & i2) : inherited(b, i1, i2), symtree(not_symmetric())
+indexed::indexed(const ex & b, const ex & i1, const ex & i2) : inherited{b, i1, i2}, symtree(not_symmetric())
 {
-       tinfo_key = &indexed::tinfo_static;
        validate();
 }
 
-indexed::indexed(const ex & b, const ex & i1, const ex & i2, const ex & i3) : inherited(b, i1, i2, i3), symtree(not_symmetric())
+indexed::indexed(const ex & b, const ex & i1, const ex & i2, const ex & i3) : inherited{b, i1, i2, i3}, symtree(not_symmetric())
 {
-       tinfo_key = &indexed::tinfo_static;
        validate();
 }
 
-indexed::indexed(const ex & b, const ex & i1, const ex & i2, const ex & i3, const ex & i4) : inherited(b, i1, i2, i3, i4), symtree(not_symmetric())
+indexed::indexed(const ex & b, const ex & i1, const ex & i2, const ex & i3, const ex & i4) : inherited{b, i1, i2, i3, i4}, symtree(not_symmetric())
 {
-       tinfo_key = &indexed::tinfo_static;
        validate();
 }
 
-indexed::indexed(const ex & b, const symmetry & symm, const ex & i1, const ex & i2) : inherited(b, i1, i2), symtree(symm)
+indexed::indexed(const ex & b, const symmetry & symm, const ex & i1, const ex & i2) : inherited{b, i1, i2}, symtree(symm)
 {
-       tinfo_key = &indexed::tinfo_static;
        validate();
 }
 
-indexed::indexed(const ex & b, const symmetry & symm, const ex & i1, const ex & i2, const ex & i3) : inherited(b, i1, i2, i3), symtree(symm)
+indexed::indexed(const ex & b, const symmetry & symm, const ex & i1, const ex & i2, const ex & i3) : inherited{b, i1, i2, i3}, symtree(symm)
 {
-       tinfo_key = &indexed::tinfo_static;
        validate();
 }
 
-indexed::indexed(const ex & b, const symmetry & symm, const ex & i1, const ex & i2, const ex & i3, const ex & i4) : inherited(b, i1, i2, i3, i4), symtree(symm)
+indexed::indexed(const ex & b, const symmetry & symm, const ex & i1, const ex & i2, const ex & i3, const ex & i4) : inherited{b, i1, i2, i3, i4}, symtree(symm)
 {
-       tinfo_key = &indexed::tinfo_static;
        validate();
 }
 
-indexed::indexed(const ex & b, const exvector & v) : inherited(b), symtree(not_symmetric())
+indexed::indexed(const ex & b, const exvector & v) : inherited{b}, symtree(not_symmetric())
 {
        seq.insert(seq.end(), v.begin(), v.end());
-       tinfo_key = &indexed::tinfo_static;
        validate();
 }
 
-indexed::indexed(const ex & b, const symmetry & symm, const exvector & v) : inherited(b), symtree(symm)
+indexed::indexed(const ex & b, const symmetry & symm, const exvector & v) : inherited{b}, symtree(symm)
 {
        seq.insert(seq.end(), v.begin(), v.end());
-       tinfo_key = &indexed::tinfo_static;
        validate();
 }
 
 indexed::indexed(const symmetry & symm, const exprseq & es) : inherited(es), symtree(symm)
 {
-       tinfo_key = &indexed::tinfo_static;
 }
 
-indexed::indexed(const symmetry & symm, const exvector & v, bool discardable) : inherited(v, discardable), symtree(symm)
+indexed::indexed(const symmetry & symm, const exvector & v) : inherited(v), symtree(symm)
 {
-       tinfo_key = &indexed::tinfo_static;
 }
 
-indexed::indexed(const symmetry & symm, std::auto_ptr<exvector> vp) : inherited(vp), symtree(symm)
+indexed::indexed(const symmetry & symm, exvector && v) : inherited(std::move(v)), symtree(symm)
 {
-       tinfo_key = &indexed::tinfo_static;
 }
 
 //////////
 // archiving
 //////////
 
-indexed::indexed(const archive_node &n, lst &sym_lst) : inherited(n, sym_lst)
+void indexed::read_archive(const archive_node &n, lst &sym_lst)
 {
+       inherited::read_archive(n, sym_lst);
        if (!n.find_ex("symmetry", symtree, sym_lst)) {
                // GiNaC versions <= 0.9.0 had an unsigned "symmetry" property
                unsigned symm = 0;
@@ -162,6 +150,7 @@ indexed::indexed(const archive_node &n, lst &sym_lst) : inherited(n, sym_lst)
                const_cast<symmetry &>(ex_to<symmetry>(symtree)).validate(seq.size() - 1);
        }
 }
+GINAC_BIND_UNARCHIVER(indexed);
 
 void indexed::archive(archive_node &n) const
 {
@@ -169,8 +158,6 @@ void indexed::archive(archive_node &n) const
        n.add_ex("symmetry", symtree);
 }
 
-DEFAULT_UNARCHIVE(indexed)
-
 //////////
 // functions overriding virtual functions from base classes
 //////////
@@ -179,7 +166,7 @@ void indexed::printindices(const print_context & c, unsigned level) const
 {
        if (seq.size() > 1) {
 
-               exvector::const_iterator it=seq.begin() + 1, itend = seq.end();
+               auto it = seq.begin() + 1, itend = seq.end();
 
                if (is_a<print_latex>(c)) {
 
@@ -256,12 +243,6 @@ bool indexed::info(unsigned inf) const
        return inherited::info(inf);
 }
 
-struct idx_is_not : public std::binary_function<ex, unsigned, bool> {
-       bool operator() (const ex & e, unsigned inf) const {
-               return !(ex_to<idx>(e).get_value().info(inf));
-       }
-};
-
 bool indexed::all_index_values_are(unsigned inf) const
 {
        // No indices? Then no property can be fulfilled
@@ -269,7 +250,8 @@ bool indexed::all_index_values_are(unsigned inf) const
                return false;
 
        // Check all indices
-       return find_if(seq.begin() + 1, seq.end(), bind2nd(idx_is_not(), inf)) == seq.end();
+       return find_if(seq.begin() + 1, seq.end(),
+                      [inf](const ex & e) { return !(ex_to<idx>(e).get_value().info(inf)); }) == seq.end();
 }
 
 int indexed::compare_same_type(const basic & other) const
@@ -298,7 +280,7 @@ ex indexed::eval(int level) const
                return f * thiscontainer(v);
        }
 
-       if(this->tinfo()==&indexed::tinfo_static && seq.size()==1)
+       if((typeid(*this) == typeid(indexed)) && seq.size()==1)
                return base;
 
        // Canonicalize indices according to the symmetry properties
@@ -306,7 +288,7 @@ ex indexed::eval(int level) const
                exvector v = seq;
                GINAC_ASSERT(is_exactly_a<symmetry>(symtree));
                int sig = canonicalize(v.begin() + 1, ex_to<symmetry>(symtree));
-               if (sig != INT_MAX) {
+               if (sig != std::numeric_limits<int>::max()) {
                        // Something has changed while sorting indices, more evaluations later
                        if (sig == 0)
                                return _ex0;
@@ -337,9 +319,9 @@ ex indexed::thiscontainer(const exvector & v) const
        return indexed(ex_to<symmetry>(symtree), v);
 }
 
-ex indexed::thiscontainer(std::auto_ptr<exvector> vp) const
+ex indexed::thiscontainer(exvector && v) const
 {
-       return indexed(ex_to<symmetry>(symtree), vp);
+       return indexed(ex_to<symmetry>(symtree), std::move(v));
 }
 
 unsigned indexed::return_type() const
@@ -390,7 +372,7 @@ ex indexed::expand(unsigned options) const
 void indexed::validate() const
 {
        GINAC_ASSERT(seq.size() > 0);
-       exvector::const_iterator it = seq.begin() + 1, itend = seq.end();
+       auto it = seq.begin() + 1, itend = seq.end();
        while (it != itend) {
                if (!is_a<idx>(*it))
                        throw(std::invalid_argument("indices of indexed object must be of type idx"));
@@ -467,7 +449,7 @@ exvector indexed::get_dummy_indices(const indexed & other) const
 
 bool indexed::has_dummy_index_for(const ex & i) const
 {
-       exvector::const_iterator it = seq.begin() + 1, itend = seq.end();
+       auto it = seq.begin() + 1, itend = seq.end();
        while (it != itend) {
                if (is_dummy_pair(*it, i))
                        return true;
@@ -545,8 +527,8 @@ exvector integral::get_free_indices() const
 template<class T> size_t number_of_type(const exvector&v)
 {
        size_t number = 0;
-       for(exvector::const_iterator i=v.begin(); i!=v.end(); ++i)
-               if(is_exactly_a<T>(*i))
+       for (auto & it : v)
+               if (is_exactly_a<T>(it))
                        ++number;
        return number;
 }
@@ -574,7 +556,7 @@ template<class T> static ex rename_dummy_indices(const ex & e, exvector & global
                // to the global set
                size_t old_global_size = global_size;
                int remaining = local_size - global_size;
-               exvector::const_iterator it = local_dummy_indices.begin(), itend = local_dummy_indices.end();
+               auto it = local_dummy_indices.begin(), itend = local_dummy_indices.end();
                while (it != itend && remaining > 0) {
                        if (is_exactly_a<T>(*it) && find_if(global_dummy_indices.begin(), global_dummy_indices.end(), bind2nd(idx_is_equal_ignore_dim(), *it)) == global_dummy_indices.end()) {
                                global_dummy_indices.push_back(*it);
@@ -639,60 +621,58 @@ bool reposition_dummy_indices(ex & e, exvector & variant_dummy_indices, exvector
 {
        bool something_changed = false;
 
-       // Canonicalize wrt indices that are dummies within e. I.e., their
-       // symbol occurs twice in an index of e. This is only done if there
-       // is a cyclic symmetry because in that case it may happen that after
-       // raising/lowering an index the indices get reshuffled by ::eval in
-       // such a way that the iterators no longer point to the right objects.
-       if (ex_to<symmetry>(ex_to<indexed>(e).get_symmetry()).has_cyclic()) {
-               // Get dummy pairs of varidxes within the indexed object in e.
-               exvector local_var_dummies;
-               local_var_dummies.reserve(e.nops()/2);
-               for (size_t i=1; i<e.nops(); ++i) {
-                       if (!is_a<varidx>(e.op(i)))
-                               continue;
-                       for (size_t j=i+1; j<e.nops(); ++j) {
-                               if (is_dummy_pair(e.op(i), e.op(j))) {
-                                       local_var_dummies.push_back(e.op(i));
-                                       for (exvector::iterator k = variant_dummy_indices.begin();
-                                                       k!=variant_dummy_indices.end(); ++k) {
-                                               if (e.op(i).op(0) == k->op(0)) {
-                                                       variant_dummy_indices.erase(k);
-                                                       break;
-                                               }
+       // Find dummy symbols that occur twice in the same indexed object.
+       exvector local_var_dummies;
+       local_var_dummies.reserve(e.nops()/2);
+       for (size_t i=1; i<e.nops(); ++i) {
+               if (!is_a<varidx>(e.op(i)))
+                       continue;
+               for (size_t j=i+1; j<e.nops(); ++j) {
+                       if (is_dummy_pair(e.op(i), e.op(j))) {
+                               local_var_dummies.push_back(e.op(i));
+                               for (auto k = variant_dummy_indices.begin(); k!=variant_dummy_indices.end(); ++k) {
+                                       if (e.op(i).op(0) == k->op(0)) {
+                                               variant_dummy_indices.erase(k);
+                                               break;
                                        }
-                                       break;
                                }
+                               break;
                        }
                }
-               // Try all posibilities of raising/lowering and keep the least one in
-               // the sense of ex_is_less.
-               ex optimal_e = e;
-               size_t numpossibs = 1 << local_var_dummies.size();
-               for (size_t i=0; i<numpossibs; ++i) {
-                       ex try_e = e;
-                       for (size_t j=0; j<local_var_dummies.size(); ++j) {
-                               exmap m;
-                               if (1<<j & i) {
-                                       ex curr_idx = local_var_dummies[j];
-                                       ex curr_toggle = ex_to<varidx>(curr_idx).toggle_variance();
-                                       m[curr_idx] = curr_toggle;
-                                       m[curr_toggle] = curr_idx;
-                               }
-                               try_e = e.subs(m, subs_options::no_pattern);
-                       }
-                       if(ex_is_less()(try_e, optimal_e))
-                       {       optimal_e = try_e;
-                               something_changed = true;
+       }
+
+       // In the case where a dummy symbol occurs twice in the same indexed object
+       // we try all possibilities of raising/lowering and keep the least one in
+       // the sense of ex_is_less.
+       ex optimal_e = e;
+       size_t numpossibs = 1 << local_var_dummies.size();
+       for (size_t i=0; i<numpossibs; ++i) {
+               ex try_e = e;
+               for (size_t j=0; j<local_var_dummies.size(); ++j) {
+                       exmap m;
+                       if (1<<j & i) {
+                               ex curr_idx = local_var_dummies[j];
+                               ex curr_toggle = ex_to<varidx>(curr_idx).toggle_variance();
+                               m[curr_idx] = curr_toggle;
+                               m[curr_toggle] = curr_idx;
                        }
+                       try_e = e.subs(m, subs_options::no_pattern);
+               }
+               if(ex_is_less()(try_e, optimal_e))
+               {       optimal_e = try_e;
+                       something_changed = true;
                }
-               e = optimal_e;
        }
+       e = optimal_e;
+
+       if (!is_a<indexed>(e))
+               return true;
+
+       exvector seq = ex_to<indexed>(e).seq;
 
        // If a dummy index is encountered for the first time in the
        // product, pull it up, otherwise, pull it down
-       exvector::const_iterator it2, it2start, it2end;
-       for (it2start = ex_to<indexed>(e).seq.begin(), it2end = ex_to<indexed>(e).seq.end(), it2 = it2start + 1; it2 != it2end; ++it2) {
+       for (auto it2 = seq.begin()+1, it2end = seq.end(); it2 != it2end; ++it2) {
                if (!is_exactly_a<varidx>(*it2))
                        continue;
 
@@ -700,14 +680,20 @@ bool reposition_dummy_indices(ex & e, exvector & variant_dummy_indices, exvector
                for (vit = variant_dummy_indices.begin(), vitend = variant_dummy_indices.end(); vit != vitend; ++vit) {
                        if (it2->op(0).is_equal(vit->op(0))) {
                                if (ex_to<varidx>(*it2).is_covariant()) {
-                                       e = e.subs(lst(
-                                               *it2 == ex_to<varidx>(*it2).toggle_variance(),
-                                               ex_to<varidx>(*it2).toggle_variance() == *it2
-                                       ), subs_options::no_pattern);
+                                       /*
+                                        * N.B. we don't want to use
+                                        *
+                                        *  e = e.subs(lst{
+                                        *  *it2 == ex_to<varidx>(*it2).toggle_variance(),
+                                        *  ex_to<varidx>(*it2).toggle_variance() == *it2
+                                        *  }, subs_options::no_pattern);
+                                        *
+                                        * since this can trigger non-trivial repositioning of indices,
+                                        * e.g. due to non-trivial symmetry properties of e, thus
+                                        * invalidating iterators
+                                        */
+                                       *it2 = ex_to<varidx>(*it2).toggle_variance();
                                        something_changed = true;
-                                       it2 = ex_to<indexed>(e).seq.begin() + (it2 - it2start);
-                                       it2start = ex_to<indexed>(e).seq.begin();
-                                       it2end = ex_to<indexed>(e).seq.end();
                                }
                                moved_indices.push_back(*vit);
                                variant_dummy_indices.erase(vit);
@@ -718,11 +704,8 @@ bool reposition_dummy_indices(ex & e, exvector & variant_dummy_indices, exvector
                for (vit = moved_indices.begin(), vitend = moved_indices.end(); vit != vitend; ++vit) {
                        if (it2->op(0).is_equal(vit->op(0))) {
                                if (ex_to<varidx>(*it2).is_contravariant()) {
-                                       e = e.subs(*it2 == ex_to<varidx>(*it2).toggle_variance(), subs_options::no_pattern);
+                                       *it2 = ex_to<varidx>(*it2).toggle_variance();
                                        something_changed = true;
-                                       it2 = ex_to<indexed>(e).seq.begin() + (it2 - it2start);
-                                       it2start = ex_to<indexed>(e).seq.begin();
-                                       it2end = ex_to<indexed>(e).seq.end();
                                }
                                goto next_index;
                        }
@@ -731,6 +714,9 @@ bool reposition_dummy_indices(ex & e, exvector & variant_dummy_indices, exvector
 next_index: ;
        }
 
+       if (something_changed)
+               e = ex_to<indexed>(e).thiscontainer(seq);
+
        return something_changed;
 }
 
@@ -778,9 +764,9 @@ static void product_to_exvector(const ex & e, exvector & v, bool & non_commutati
 template<class T> ex idx_symmetrization(const ex& r,const exvector& local_dummy_indices)
 {      exvector dummy_syms;
        dummy_syms.reserve(r.nops());
-       for (exvector::const_iterator it = local_dummy_indices.begin(); it != local_dummy_indices.end(); ++it)
-                       if(is_exactly_a<T>(*it))
-                               dummy_syms.push_back(it->op(0));
+       for (auto & it : local_dummy_indices)
+                       if(is_exactly_a<T>(it))
+                               dummy_syms.push_back(it.op(0));
        if(dummy_syms.size() < 2)
                return r;
        ex q=symmetrize(r, dummy_syms);
@@ -804,6 +790,7 @@ ex simplify_indexed_product(const ex & e, exvector & free_indices, exvector & du
 
        // Perform contractions
        bool something_changed = false;
+       bool has_nonsymmetric = false;
        GINAC_ASSERT(v.size() > 1);
        exvector::iterator it1, itend = v.end(), next_to_last = itend - 1;
        for (it1 = v.begin(); it1 != next_to_last; it1++) {
@@ -813,6 +800,7 @@ try_again:
                        continue;
 
                bool first_noncommutative = (it1->return_type() != return_types::commutative);
+               bool first_nonsymmetric = ex_to<symmetry>(ex_to<indexed>(*it1).get_symmetry()).has_nonsymmetric();
 
                // Indexed factor found, get free indices and look for contraction
                // candidates
@@ -879,8 +867,16 @@ contraction_done:
                                        // Non-commutative products are always re-expanded to give
                                        // eval_ncmul() the chance to re-order and canonicalize
                                        // the product
-                                       ex r = (non_commutative ? ex(ncmul(v, true)) : ex(mul(v)));
-                                       return simplify_indexed(r, free_indices, dummy_indices, sp);
+                                       bool is_a_product = (is_exactly_a<mul>(*it1) || is_exactly_a<ncmul>(*it1)) &&
+                                                           (is_exactly_a<mul>(*it2) || is_exactly_a<ncmul>(*it2));
+                                       ex r = (non_commutative ? ex(ncmul(std::move(v))) : ex(mul(std::move(v))));
+
+                                       // If new expression is a product we can call this function again,
+                                       // otherwise we need to pass argument to simplify_indexed() to be expanded
+                                       if (is_a_product)
+                                               return simplify_indexed_product(r, free_indices, dummy_indices, sp);
+                                       else
+                                               return simplify_indexed(r, free_indices, dummy_indices, sp);
                                }
 
                                // Both objects may have new indices now or they might
@@ -889,6 +885,11 @@ contraction_done:
                                something_changed = true;
                                goto try_again;
                        }
+                       else if (!has_nonsymmetric &&
+                                       (first_nonsymmetric ||
+                                        ex_to<symmetry>(ex_to<indexed>(*it2).get_symmetry()).has_nonsymmetric())) {
+                               has_nonsymmetric = true;
+                       }
                }
        }
 
@@ -935,27 +936,29 @@ contraction_done:
 
        ex r;
        if (something_changed)
-               r = non_commutative ? ex(ncmul(v, true)) : ex(mul(v));
+               r = non_commutative ? ex(ncmul(std::move(v))) : ex(mul(std::move(v)));
        else
                r = e;
 
        // The result should be symmetric with respect to exchange of dummy
        // indices, so if the symmetrization vanishes, the whole expression is
        // zero. This detects things like eps.i.j.k * p.j * p.k = 0.
-       ex q = idx_symmetrization<idx>(r, local_dummy_indices);
-       if (q.is_zero()) {
-               free_indices.clear();
-               return _ex0;
-       }
-       q = idx_symmetrization<varidx>(q, local_dummy_indices);
-       if (q.is_zero()) {
-               free_indices.clear();
-               return _ex0;
-       }
-       q = idx_symmetrization<spinidx>(q, local_dummy_indices);
-       if (q.is_zero()) {
-               free_indices.clear();
-               return _ex0;
+       if (has_nonsymmetric) {
+               ex q = idx_symmetrization<idx>(r, local_dummy_indices);
+               if (q.is_zero()) {
+                       free_indices.clear();
+                       return _ex0;
+               }
+               q = idx_symmetrization<varidx>(q, local_dummy_indices);
+               if (q.is_zero()) {
+                       free_indices.clear();
+                       return _ex0;
+               }
+               q = idx_symmetrization<spinidx>(q, local_dummy_indices);
+               if (q.is_zero()) {
+                       free_indices.clear();
+                       return _ex0;
+               }
        }
 
        // Dummy index renaming
@@ -978,7 +981,7 @@ public:
        terminfo(const ex & orig_, const ex & symm_) : orig(orig_), symm(symm_) {}
 
        ex orig; /**< original term */
-       ex symm; /**< symmtrized term */
+       ex symm; /**< symmetrized term */
 };
 
 class terminfo_is_less {
@@ -1121,9 +1124,9 @@ ex simplify_indexed(const ex & e, exvector & free_indices, exvector & dummy_indi
                        const ex & term = sum.op(i);
                        exvector dummy_indices_of_term;
                        dummy_indices_of_term.reserve(dummy_indices.size());
-                       for(exvector::iterator i=dummy_indices.begin(); i!=dummy_indices.end(); ++i)
-                               if(hasindex(term,i->op(0)))
-                                       dummy_indices_of_term.push_back(*i);
+                       for (auto & i : dummy_indices)
+                               if (hasindex(term,i.op(0)))
+                                       dummy_indices_of_term.push_back(i);
                        ex term_symm = idx_symmetrization<idx>(term, dummy_indices_of_term);
                        term_symm = idx_symmetrization<varidx>(term_symm, dummy_indices_of_term);
                        term_symm = idx_symmetrization<spinidx>(term_symm, dummy_indices_of_term);
@@ -1139,7 +1142,7 @@ ex simplify_indexed(const ex & e, exvector & free_indices, exvector & dummy_indi
                std::vector<terminfo> terms_pass2;
                for (std::vector<terminfo>::const_iterator i=terms.begin(); i!=terms.end(); ) {
                        size_t num = 1;
-                       std::vector<terminfo>::const_iterator j = i + 1;
+                       auto j = i + 1;
                        while (j != terms.end() && j->symm == i->symm) {
                                num++;
                                j++;
@@ -1154,13 +1157,13 @@ ex simplify_indexed(const ex & e, exvector & free_indices, exvector & dummy_indi
 
                // Chop the symmetrized terms into subterms
                std::vector<symminfo> sy;
-               for (std::vector<terminfo>::const_iterator i=terms_pass2.begin(); i!=terms_pass2.end(); ++i) {
-                       if (is_exactly_a<add>(i->symm)) {
-                               size_t num = i->symm.nops();
+               for (auto & i : terms_pass2) {
+                       if (is_exactly_a<add>(i.symm)) {
+                               size_t num = i.symm.nops();
                                for (size_t j=0; j<num; j++)
-                                       sy.push_back(symminfo(i->symm.op(j), i->orig, num));
+                                       sy.push_back(symminfo(i.symm.op(j), i.orig, num));
                        } else
-                               sy.push_back(symminfo(i->symm, i->orig, 1));
+                               sy.push_back(symminfo(i.symm, i.orig, 1));
                }
 
                // Sort by symmetrized subterms
@@ -1169,10 +1172,10 @@ ex simplify_indexed(const ex & e, exvector & free_indices, exvector & dummy_indi
                // Combine equal symmetrized subterms
                std::vector<symminfo> sy_pass2;
                exvector result;
-               for (std::vector<symminfo>::const_iterator i=sy.begin(); i!=sy.end(); ) {
+               for (auto i=sy.begin(); i!=sy.end(); ) {
 
                        // Combine equal terms
-                       std::vector<symminfo>::const_iterator j = i + 1;
+                       auto j = i + 1;
                        if (j != sy.end() && j->symmterm == i->symmterm) {
 
                                // More than one term, collect the coefficients
@@ -1205,7 +1208,7 @@ ex simplify_indexed(const ex & e, exvector & free_indices, exvector & dummy_indi
 
                                // How many symmetrized terms of this original term are left?
                                size_t num = 1;
-                               std::vector<symminfo>::const_iterator j = i + 1;
+                               auto j = i + 1;
                                while (j != sy_pass2.end() && j->orig == i->orig) {
                                        num++;
                                        j++;
@@ -1229,7 +1232,7 @@ ex simplify_indexed(const ex & e, exvector & free_indices, exvector & dummy_indi
                }
 
                // Add all resulting terms
-               ex sum_symm = (new add(result))->setflag(status_flags::dynallocated);
+               ex sum_symm = dynallocate<add>(result);
                if (sum_symm.is_zero())
                        free_indices.clear();
                return sum_symm;
@@ -1357,9 +1360,9 @@ void scalar_products::add(const ex & v1, const ex & v2, const ex & dim, const ex
 void scalar_products::add_vectors(const lst & l, const ex & dim)
 {
        // Add all possible pairs of products
-       for (lst::const_iterator it1 = l.begin(); it1 != l.end(); ++it1)
-               for (lst::const_iterator it2 = l.begin(); it2 != l.end(); ++it2)
-                       add(*it1, *it2, *it1 * *it2);
+       for (auto & it1 : l)
+               for (auto & it2 : l)
+                       add(it1, it2, it1 * it2);
 }
 
 void scalar_products::clear()
@@ -1382,13 +1385,11 @@ ex scalar_products::evaluate(const ex & v1, const ex & v2, const ex & dim) const
 void scalar_products::debugprint() const
 {
        std::cerr << "map size=" << spm.size() << std::endl;
-       spmap::const_iterator i = spm.begin(), end = spm.end();
-       while (i != end) {
-               const spmapkey & k = i->first;
+       for (auto & it : spm) {
+               const spmapkey & k = it.first;
                std::cerr << "item key=";
                k.debugprint();
-               std::cerr << ", value=" << i->second << std::endl;
-               ++i;
+               std::cerr << ", value=" << it.second << std::endl;
        }
 }
 
@@ -1402,7 +1403,7 @@ exvector get_all_dummy_indices_safely(const ex & e)
        else if (is_a<mul>(e) || is_a<ncmul>(e)) {
                exvector dummies;
                exvector free_indices;
-               for (int i=0; i<e.nops(); ++i) {
+               for (std::size_t i = 0; i < e.nops(); ++i) {
                        exvector dummies_of_factor = get_all_dummy_indices_safely(e.op(i));
                        dummies.insert(dummies.end(), dummies_of_factor.begin(),
                                dummies_of_factor.end());
@@ -1418,7 +1419,7 @@ exvector get_all_dummy_indices_safely(const ex & e)
        }
        else if(is_a<add>(e)) {
                exvector result;
-               for(int i=0; i<e.nops(); ++i) {
+               for(std::size_t i = 0; i < e.nops(); ++i) {
                        exvector dummies_of_term = get_all_dummy_indices_safely(e.op(i));
                        sort(dummies_of_term.begin(), dummies_of_term.end());
                        exvector new_vec;
@@ -1438,13 +1439,13 @@ exvector get_all_dummy_indices(const ex & e)
        exvector p;
        bool nc;
        product_to_exvector(e, p, nc);
-       exvector::const_iterator ip = p.begin(), ipend = p.end();
+       auto ip = p.begin(), ipend = p.end();
        exvector v, v1;
        while (ip != ipend) {
                if (is_a<indexed>(*ip)) {
                        v1 = ex_to<indexed>(*ip).get_dummy_indices();
                        v.insert(v.end(), v1.begin(), v1.end());
-                       exvector::const_iterator ip1 = ip+1;
+                       auto ip1 = ip + 1;
                        while (ip1 != ipend) {
                                if (is_a<indexed>(*ip1)) {
                                        v1 = ex_to<indexed>(*ip).get_dummy_indices(ex_to<indexed>(*ip1));
@@ -1463,27 +1464,24 @@ lst rename_dummy_indices_uniquely(const exvector & va, const exvector & vb)
        exvector common_indices;
        set_intersection(va.begin(), va.end(), vb.begin(), vb.end(), std::back_insert_iterator<exvector>(common_indices), ex_is_less());
        if (common_indices.empty()) {
-               return lst(lst(), lst());
+               return lst{lst{}, lst{}};
        } else {
                exvector new_indices, old_indices;
                old_indices.reserve(2*common_indices.size());
                new_indices.reserve(2*common_indices.size());
                exvector::const_iterator ip = common_indices.begin(), ipend = common_indices.end();
                while (ip != ipend) {
-                       ex newsym=(new symbol)->setflag(status_flags::dynallocated);
+                       ex newsym = dynallocate<symbol>();
                        ex newidx;
                        if(is_exactly_a<spinidx>(*ip))
-                               newidx = (new spinidx(newsym, ex_to<spinidx>(*ip).get_dim(),
-                                               ex_to<spinidx>(*ip).is_covariant(),
-                                               ex_to<spinidx>(*ip).is_dotted()))
-                                       -> setflag(status_flags::dynallocated);
+                               newidx = dynallocate<spinidx>(newsym, ex_to<spinidx>(*ip).get_dim(),
+                                                             ex_to<spinidx>(*ip).is_covariant(),
+                                                             ex_to<spinidx>(*ip).is_dotted());
                        else if (is_exactly_a<varidx>(*ip))
-                               newidx = (new varidx(newsym, ex_to<varidx>(*ip).get_dim(),
-                                               ex_to<varidx>(*ip).is_covariant()))
-                                       -> setflag(status_flags::dynallocated);
+                               newidx = dynallocate<varidx>(newsym, ex_to<varidx>(*ip).get_dim(),
+                                                            ex_to<varidx>(*ip).is_covariant());
                        else
-                               newidx = (new idx(newsym, ex_to<idx>(*ip).get_dim()))
-                                       -> setflag(status_flags::dynallocated);
+                               newidx = dynallocate<idx>(newsym, ex_to<idx>(*ip).get_dim());
                        old_indices.push_back(*ip);
                        new_indices.push_back(newidx);
                        if(is_a<varidx>(*ip)) {
@@ -1492,7 +1490,7 @@ lst rename_dummy_indices_uniquely(const exvector & va, const exvector & vb)
                        }
                        ++ip;
                }
-               return lst(lst(old_indices.begin(), old_indices.end()), lst(new_indices.begin(), new_indices.end()));
+               return lst{lst(old_indices.begin(), old_indices.end()), lst(new_indices.begin(), new_indices.end())};
        }
 }
 
@@ -1527,15 +1525,12 @@ ex rename_dummy_indices_uniquely(exvector & va, const ex & b, bool modify_va)
                        lst indices_subs = rename_dummy_indices_uniquely(va, vb);
                        if (indices_subs.op(0).nops() > 0) {
                                if (modify_va) {
-                                       for (lst::const_iterator i = ex_to<lst>(indices_subs.op(1)).begin(); i != ex_to<lst>(indices_subs.op(1)).end(); ++i)
-                                               va.push_back(*i);
+                                       for (auto & i : ex_to<lst>(indices_subs.op(1)))
+                                               va.push_back(i);
                                        exvector uncommon_indices;
                                        set_difference(vb.begin(), vb.end(), indices_subs.op(0).begin(), indices_subs.op(0).end(), std::back_insert_iterator<exvector>(uncommon_indices), ex_is_less());
-                                       exvector::const_iterator ip = uncommon_indices.begin(), ipend = uncommon_indices.end();
-                                       while (ip != ipend) {
-                                               va.push_back(*ip);
-                                               ++ip;
-                                       }
+                                       for (auto & ip : uncommon_indices)
+                                               va.push_back(ip);
                                        sort(va.begin(), va.end(), ex_is_less());
                                }
                                return b.subs(ex_to<lst>(indices_subs.op(0)), ex_to<lst>(indices_subs.op(1)), subs_options::no_pattern|subs_options::no_index_renaming);
@@ -1558,18 +1553,17 @@ ex expand_dummy_sum(const ex & e, bool subs_idx)
                else
                        v = get_all_dummy_indices(e_expanded);
                ex result = e_expanded;
-               for(exvector::const_iterator it=v.begin(); it!=v.end(); ++it) {
-                       ex nu = *it;
+               for (const auto & nu : v) {
                        if (ex_to<idx>(nu).get_dim().info(info_flags::nonnegint)) {
                                int idim = ex_to<numeric>(ex_to<idx>(nu).get_dim()).to_int();
                                ex en = 0;
                                for (int i=0; i < idim; i++) {
                                        if (subs_idx && is_a<varidx>(nu)) {
                                                ex other = ex_to<varidx>(nu).toggle_variance();
-                                               en += result.subs(lst(
+                                               en += result.subs(lst{
                                                        nu == idx(i, idim),
                                                        other == idx(i, idim)
-                                               ));
+                                               });
                                        } else {
                                                en += result.subs( nu.op(0) == i );
                                        }