Index renaming issues, sped up simplify_indexed, used defined NC-objects

can carry indices.

diff --git a/doc/tutorial/ginac.texi b/doc/tutorial/ginac.texi
index 06ab56fe55ae1c9baed9498d7d2bbcf48d35de33..3cd6c6179478ce9d05b81261eb9317e53b639744 100644 (file)
--- a/doc/tutorial/ginac.texi
+++ b/doc/tutorial/ginac.texi
@@ -4278,13 +4278,17 @@ contain the same number of elements). Using this form, you would write
 @end example
 
 The optional last argument to @code{subs()} is a combination of
 @end example
 
 The optional last argument to @code{subs()} is a combination of

-@code{subs_options} flags. There are two options available:
+@code{subs_options} flags. There are three options available:

 @code{subs_options::no_pattern} disables pattern matching, which makes
 large @code{subs()} operations significantly faster if you are not using
 patterns. The second option, @code{subs_options::algebraic} enables
 algebraic substitutions in products and powers.
 @ref{Pattern Matching and Advanced Substitutions}, for more information
 @code{subs_options::no_pattern} disables pattern matching, which makes
 large @code{subs()} operations significantly faster if you are not using
 patterns. The second option, @code{subs_options::algebraic} enables
 algebraic substitutions in products and powers.
 @ref{Pattern Matching and Advanced Substitutions}, for more information

-about patterns and algebraic substitutions.
+about patterns and algebraic substitutions. The third option,
+@code{subs_options::no_index_renaming} disables the feature that dummy
+indices are renamed if the subsitution could give a result in which a
+dummy index occurs more than two times. This is sometimes necessary if
+you want to use @code{subs()} to rename your dummy indices.

 
 @code{subs()} performs syntactic substitution of any complete algebraic
 object; it does not try to match sub-expressions as is demonstrated by the
 
 @code{subs()} performs syntactic substitution of any complete algebraic
 object; it does not try to match sub-expressions as is demonstrated by the

diff --git a/ginac/color.h b/ginac/color.h
index a585f2e575f3b08af460b80f70329a3f1826c8dc..9d9ba26e88e1e540b914f24185e2814e27493bd7 100644 (file)
--- a/ginac/color.h
+++ b/ginac/color.h
@@ -109,6 +109,7 @@ public:
 
        // non-virtual functions in this class
 protected:
 
        // non-virtual functions in this class
 protected:

+       unsigned return_type() const { return return_types::commutative; }

        void do_print(const print_context & c, unsigned level) const;
        void do_print_latex(const print_latex & c, unsigned level) const;
 };
        void do_print(const print_context & c, unsigned level) const;
        void do_print_latex(const print_latex & c, unsigned level) const;
 };


@@ -125,6 +126,7 @@ public:
 
        // non-virtual functions in this class
 protected:
 
        // non-virtual functions in this class
 protected:

+       unsigned return_type() const { return return_types::commutative; }

        void do_print(const print_context & c, unsigned level) const;
        void do_print_latex(const print_latex & c, unsigned level) const;
 };
        void do_print(const print_context & c, unsigned level) const;
        void do_print_latex(const print_latex & c, unsigned level) const;
 };

diff --git a/ginac/expairseq.cpp b/ginac/expairseq.cpp
index 7555bd185d422c011beb9a86f5de8e6ea5b61666..9f580d9e768ff4ac58c5fcaaa0865124ad7069fd 100644 (file)
--- a/ginac/expairseq.cpp
+++ b/ginac/expairseq.cpp
@@ -34,6 +34,7 @@
 #include "archive.h"
 #include "operators.h"
 #include "utils.h"
 #include "archive.h"
 #include "operators.h"
 #include "utils.h"

+#include "indexed.h"

 
 #if EXPAIRSEQ_USE_HASHTAB
 #include <cmath>
 
 #if EXPAIRSEQ_USE_HASHTAB
 #include <cmath>


@@ -757,8 +758,15 @@ void expairseq::construct_from_2_ex(const ex &lh, const ex &rh)
                                construct_from_2_ex_via_exvector(lh,rh);
                        } else {
 #endif // EXPAIRSEQ_USE_HASHTAB
                                construct_from_2_ex_via_exvector(lh,rh);
                        } else {
 #endif // EXPAIRSEQ_USE_HASHTAB

-                               construct_from_2_expairseq(ex_to<expairseq>(lh),
-                                                          ex_to<expairseq>(rh));
+                               if(is_a<mul>(lh))
+                               {       
+                                       ex newrh=rename_dummy_indices_uniquely(lh, rh);
+                                       construct_from_2_expairseq(ex_to<expairseq>(lh),
+                                                                  ex_to<expairseq>(newrh));
+                               }
+                               else
+                                       construct_from_2_expairseq(ex_to<expairseq>(lh),
+                                                                  ex_to<expairseq>(rh));

 #if EXPAIRSEQ_USE_HASHTAB
                        }
 #endif // EXPAIRSEQ_USE_HASHTAB
 #if EXPAIRSEQ_USE_HASHTAB
                        }
 #endif // EXPAIRSEQ_USE_HASHTAB


@@ -1008,13 +1016,27 @@ void expairseq::make_flat(const exvector &v)
        seq.reserve(v.size()+noperands-nexpairseqs);
        
        // copy elements and split off numerical part
        seq.reserve(v.size()+noperands-nexpairseqs);
        
        // copy elements and split off numerical part

+       exvector dummy_indices;

        cit = v.begin();
        while (cit!=v.end()) {
                if (ex_to<basic>(*cit).tinfo()==this->tinfo()) {
        cit = v.begin();
        while (cit!=v.end()) {
                if (ex_to<basic>(*cit).tinfo()==this->tinfo()) {

-                       const expairseq &subseqref = ex_to<expairseq>(*cit);
-                       combine_overall_coeff(subseqref.overall_coeff);
-                       epvector::const_iterator cit_s = subseqref.seq.begin();
-                       while (cit_s!=subseqref.seq.end()) {
+                       const expairseq *subseqref;
+                       ex newfactor;
+                       if(is_a<mul>(*cit))
+                       {
+                               exvector dummies_of_factor = get_all_dummy_indices(*cit);
+                               sort(dummies_of_factor.begin(), dummies_of_factor.end(), ex_is_less());
+                               newfactor = rename_dummy_indices_uniquely(dummy_indices, dummies_of_factor, *cit);
+                               subseqref = &(ex_to<expairseq>(newfactor));
+                               exvector new_dummy_indices;
+                               set_union(dummy_indices.begin(), dummy_indices.end(), dummies_of_factor.begin(), dummies_of_factor.end(), std::back_insert_iterator<exvector>(new_dummy_indices), ex_is_less());
+                               dummy_indices.swap(new_dummy_indices);
+                       }
+                       else
+                               subseqref = &ex_to<expairseq>(*cit);
+                       combine_overall_coeff(subseqref->overall_coeff);
+                       epvector::const_iterator cit_s = subseqref->seq.begin();
+                       while (cit_s!=subseqref->seq.end()) {

                                seq.push_back(*cit_s);
                                ++cit_s;
                        }
                                seq.push_back(*cit_s);
                                ++cit_s;
                        }


@@ -1579,6 +1601,67 @@ std::auto_ptr<epvector> expairseq::evalchildren(int level) const
        return std::auto_ptr<epvector>(0); // signalling nothing has changed
 }
 
        return std::auto_ptr<epvector>(0); // signalling nothing has changed
 }
 

+class safe_inserter
+{
+       public:
+               safe_inserter(const ex&, const bool disable_renaming=false);
+               std::auto_ptr<epvector> getseq(){return epv;}
+               void insert_old_pair(const expair &p)
+               {
+                       epv->push_back(p);
+               }
+               void insert_new_pair(const expair &p, const ex &orig_ex);
+       private:
+               std::auto_ptr<epvector> epv;
+               bool dodummies;
+               exvector dummy_indices;
+               void update_dummy_indices(const exvector&);
+};
+
+safe_inserter::safe_inserter(const ex&e, const bool disable_renaming)
+               :epv(new epvector)
+{
+       epv->reserve(e.nops());
+       dodummies=is_a<mul>(e);
+       if(disable_renaming)
+               dodummies=false;
+       if(dodummies) {
+               dummy_indices = get_all_dummy_indices(e);
+               sort(dummy_indices.begin(), dummy_indices.end(), ex_is_less());
+       }
+}
+
+void safe_inserter::update_dummy_indices(const exvector &v)
+{
+       exvector new_dummy_indices;
+       set_union(dummy_indices.begin(), dummy_indices.end(), v.begin(), v.end(),
+               std::back_insert_iterator<exvector>(new_dummy_indices), ex_is_less());
+       dummy_indices.swap(new_dummy_indices);
+}
+
+void safe_inserter::insert_new_pair(const expair &p, const ex &orig_ex)
+{
+       if(!dodummies) {
+               epv->push_back(p);
+               return;
+       }
+       exvector dummies_of_factor = get_all_dummy_indices(p.rest);
+       if(dummies_of_factor.size() == 0) {
+               epv->push_back(p);
+               return;
+       }
+       sort(dummies_of_factor.begin(), dummies_of_factor.end(), ex_is_less());
+       exvector dummies_of_orig_ex = get_all_dummy_indices(orig_ex);
+       sort(dummies_of_orig_ex.begin(), dummies_of_orig_ex.end(), ex_is_less());
+       exvector new_dummy_indices;
+       new_dummy_indices.reserve(dummy_indices.size());
+       set_difference(dummy_indices.begin(), dummy_indices.end(), dummies_of_orig_ex.begin(), dummies_of_orig_ex.end(),
+               std::back_insert_iterator<exvector>(new_dummy_indices), ex_is_less());
+       dummy_indices.swap(new_dummy_indices);
+       ex newfactor = rename_dummy_indices_uniquely(dummy_indices, dummies_of_factor, p.rest);
+       update_dummy_indices(dummies_of_factor);
+       epv -> push_back(expair(newfactor, p.coeff));
+}

 
 /** Member-wise substitute in this sequence.
  *
 
 /** Member-wise substitute in this sequence.
  *


@@ -1614,22 +1697,27 @@ std::auto_ptr<epvector> expairseq::subschildren(const exmap & m, unsigned option
                        if (!are_ex_trivially_equal(orig_ex, subsed_ex)) {
 
                                // Something changed, copy seq, subs and return it
                        if (!are_ex_trivially_equal(orig_ex, subsed_ex)) {
 
                                // Something changed, copy seq, subs and return it

-                               std::auto_ptr<epvector> s(new epvector);
-                               s->reserve(seq.size());
+                               safe_inserter s(*this, options & subs_options::no_index_renaming);

 
                                // Copy parts of seq which are known not to have changed
 
                                // Copy parts of seq which are known not to have changed

-                               s->insert(s->begin(), seq.begin(), cit);
+                               for(epvector::const_iterator i=seq.begin(); i!=cit; ++i)
+                                       s.insert_old_pair(*i);

 
                                // Copy first changed element
 
                                // Copy first changed element

-                               s->push_back(split_ex_to_pair(subsed_ex));
+                               s.insert_new_pair(split_ex_to_pair(subsed_ex), orig_ex);

                                ++cit;
 
                                // Copy rest
                                while (cit != last) {
                                ++cit;
 
                                // Copy rest
                                while (cit != last) {

-                                       s->push_back(split_ex_to_pair(recombine_pair_to_ex(*cit).subs(m, options)));
+                                       ex orig_ex = recombine_pair_to_ex(*cit);
+                                       ex subsed_ex = orig_ex.subs(m, options);
+                                       if(are_ex_trivially_equal(orig_ex, subsed_ex))
+                                               s.insert_old_pair(*cit);
+                                       else
+                                               s.insert_new_pair(split_ex_to_pair(subsed_ex), orig_ex);

                                        ++cit;
                                }
                                        ++cit;
                                }

-                               return s;
+                               return s.getseq();

                        }
 
                        ++cit;
                        }
 
                        ++cit;


@@ -1645,23 +1733,27 @@ std::auto_ptr<epvector> expairseq::subschildren(const exmap & m, unsigned option
                        if (!are_ex_trivially_equal(cit->rest, subsed_ex)) {
                        
                                // Something changed, copy seq, subs and return it
                        if (!are_ex_trivially_equal(cit->rest, subsed_ex)) {
                        
                                // Something changed, copy seq, subs and return it

-                               std::auto_ptr<epvector> s(new epvector);
-                               s->reserve(seq.size());
+                               safe_inserter s(*this, options & subs_options::no_index_renaming);

 
                                // Copy parts of seq which are known not to have changed
 
                                // Copy parts of seq which are known not to have changed

-                               s->insert(s->begin(), seq.begin(), cit);
+                               for(epvector::const_iterator i=seq.begin(); i!=cit; ++i)
+                                       s.insert_old_pair(*i);

                        
                                // Copy first changed element
                        
                                // Copy first changed element

-                               s->push_back(combine_ex_with_coeff_to_pair(subsed_ex, cit->coeff));
+                               s.insert_new_pair(combine_ex_with_coeff_to_pair(subsed_ex, cit->coeff), cit->rest);

                                ++cit;
 
                                // Copy rest
                                while (cit != last) {
                                ++cit;
 
                                // Copy rest
                                while (cit != last) {

-                                       s->push_back(combine_ex_with_coeff_to_pair(cit->rest.subs(m, options),
-                                                                                  cit->coeff));
+                                       const ex &orig_ex = cit->rest;
+                                       const ex &subsed_ex = cit->rest.subs(m, options);
+                                       if(are_ex_trivially_equal(orig_ex, subsed_ex))
+                                               s.insert_old_pair(*cit);
+                                       else
+                                               s.insert_new_pair(combine_ex_with_coeff_to_pair(subsed_ex, cit->coeff), orig_ex);

                                        ++cit;
                                }
                                        ++cit;
                                }

-                               return s;
+                               return s.getseq();

                        }
 
                        ++cit;
                        }
 
                        ++cit;

diff --git a/ginac/flags.h b/ginac/flags.h
index 90f58ad2e734669dd383b7a2c7d77a2b3c3942dd..7c51b75d2ecd7d29cb29dcb8c3b3fe35240e9208 100644 (file)
--- a/ginac/flags.h
+++ b/ginac/flags.h
@@ -43,7 +43,8 @@ public:
                algebraic = 0x0002,              ///< enable algebraic substitutions
                subs_algebraic = 0x0002,  // for backwards compatibility
                pattern_is_product = 0x0004,     ///< used internally by expairseq::subschildren()
                algebraic = 0x0002,              ///< enable algebraic substitutions
                subs_algebraic = 0x0002,  // for backwards compatibility
                pattern_is_product = 0x0004,     ///< used internally by expairseq::subschildren()

-               pattern_is_not_product = 0x0008  ///< used internally by expairseq::subschildren()
+               pattern_is_not_product = 0x0008, ///< used internally by expairseq::subschildren()
+               no_index_renaming = 0x0010

        };
 };
 
        };
 };
 

diff --git a/ginac/indexed.cpp b/ginac/indexed.cpp
index ef904f7222bfc798420dee5cdb6bed6373426a26..64eba6d08cf1cccf6dbf455cbcce4c01356dd3e7 100644 (file)
--- a/ginac/indexed.cpp
+++ b/ginac/indexed.cpp
@@ -297,6 +297,9 @@ ex indexed::eval(int level) const
                return f * thiscontainer(v);
        }
 
                return f * thiscontainer(v);
        }
 

+       if(this->tinfo()==TINFO_indexed && seq.size()==1)
+               return base;
+

        // Canonicalize indices according to the symmetry properties
        if (seq.size() > 2) {
                exvector v = seq;
        // Canonicalize indices according to the symmetry properties
        if (seq.size() > 2) {
                exvector v = seq;


@@ -324,6 +327,14 @@ ex indexed::thiscontainer(std::auto_ptr<exvector> vp) const
        return indexed(ex_to<symmetry>(symtree), vp);
 }
 
        return indexed(ex_to<symmetry>(symtree), vp);
 }
 

+unsigned indexed::return_type() const
+{
+       if(is_a<matrix>(op(0)))
+               return return_types::commutative;
+       else
+               return op(0).return_type();
+}
+

 ex indexed::expand(unsigned options) const
 {
        GINAC_ASSERT(seq.size() > 0);
 ex indexed::expand(unsigned options) const
 {
        GINAC_ASSERT(seq.size() > 0);


@@ -532,6 +543,15 @@ exvector integral::get_free_indices() const
        return f.get_free_indices();
 }
 
        return f.get_free_indices();
 }
 

+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))
+                       ++number;
+       return number;
+}
+

 /** Rename dummy indices in an expression.
  *
  *  @param e Expression to work on
 /** Rename dummy indices in an expression.
  *
  *  @param e Expression to work on


@@ -540,10 +560,10 @@ exvector integral::get_free_indices() const
  *  @param global_dummy_indices The set of dummy indices that have appeared
  *    before and which we would like to use in "e", too. This gets updated
  *    by the function */
  *  @param global_dummy_indices The set of dummy indices that have appeared
  *    before and which we would like to use in "e", too. This gets updated
  *    by the function */

-static ex rename_dummy_indices(const ex & e, exvector & global_dummy_indices, exvector & local_dummy_indices)
+template<class T> static ex rename_dummy_indices(const ex & e, exvector & global_dummy_indices, exvector & local_dummy_indices)

 {
 {

-       size_t global_size = global_dummy_indices.size(),
-              local_size = local_dummy_indices.size();
+       size_t global_size = number_of_type<T>(global_dummy_indices),
+              local_size = number_of_type<T>(local_dummy_indices);

 
        // Any local dummy indices at all?
        if (local_size == 0)
 
        // Any local dummy indices at all?
        if (local_size == 0)


@@ -557,7 +577,7 @@ static ex rename_dummy_indices(const ex & e, exvector & global_dummy_indices, ex
                int remaining = local_size - global_size;
                exvector::const_iterator it = local_dummy_indices.begin(), itend = local_dummy_indices.end();
                while (it != itend && remaining > 0) {
                int remaining = local_size - global_size;
                exvector::const_iterator it = local_dummy_indices.begin(), itend = local_dummy_indices.end();
                while (it != itend && remaining > 0) {

-                       if (find_if(global_dummy_indices.begin(), global_dummy_indices.end(), bind2nd(op0_is_equal(), *it)) == global_dummy_indices.end()) {
+                       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);
                                global_size++;
                                remaining--;
                                global_dummy_indices.push_back(*it);
                                global_size++;
                                remaining--;


@@ -575,11 +595,13 @@ static ex rename_dummy_indices(const ex & e, exvector & global_dummy_indices, ex
        exvector local_syms, global_syms;
        local_syms.reserve(local_size);
        global_syms.reserve(local_size);
        exvector local_syms, global_syms;
        local_syms.reserve(local_size);
        global_syms.reserve(local_size);

-       for (size_t i=0; i<local_size; i++)
-               local_syms.push_back(local_dummy_indices[i].op(0));
+       for (size_t i=0; local_syms.size()!=local_size; i++)
+               if(is_exactly_a<T>(local_dummy_indices[i]))
+                       local_syms.push_back(local_dummy_indices[i].op(0));

        shaker_sort(local_syms.begin(), local_syms.end(), ex_is_less(), ex_swap());
        shaker_sort(local_syms.begin(), local_syms.end(), ex_is_less(), ex_swap());

-       for (size_t i=0; i<local_size; i++) // don't use more global symbols than necessary
-               global_syms.push_back(global_dummy_indices[i].op(0));
+       for (size_t i=0; global_syms.size()!=local_size; i++) // don't use more global symbols than necessary
+               if(is_exactly_a<T>(global_dummy_indices[i]))
+                       global_syms.push_back(global_dummy_indices[i].op(0));

        shaker_sort(global_syms.begin(), global_syms.end(), ex_is_less(), ex_swap());
 
        // Remove common indices
        shaker_sort(global_syms.begin(), global_syms.end(), ex_is_less(), ex_swap());
 
        // Remove common indices


@@ -704,6 +726,18 @@ 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));
+       if(dummy_syms.size() < 2)
+               return r;
+       ex q=symmetrize(r, dummy_syms);
+       return q;
+}
+

 /** Simplify product of indexed expressions (commutative, noncommutative and
  *  simple squares), return list of free indices. */
 ex simplify_indexed_product(const ex & e, exvector & free_indices, exvector & dummy_indices, const scalar_products & sp)
 /** Simplify product of indexed expressions (commutative, noncommutative and
  *  simple squares), return list of free indices. */
 ex simplify_indexed_product(const ex & e, exvector & free_indices, exvector & dummy_indices, const scalar_products & sp)


@@ -864,19 +898,26 @@ contraction_done:
        // 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.
        // 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.

-       if (local_dummy_indices.size() >= 2) {
-               exvector dummy_syms;
-               dummy_syms.reserve(local_dummy_indices.size());
-               for (exvector::const_iterator it = local_dummy_indices.begin(); it != local_dummy_indices.end(); ++it)
-                       dummy_syms.push_back(it->op(0));
-               if (symmetrize(r, dummy_syms).is_zero()) {
-                       free_indices.clear();
-                       return _ex0;
-               }
+       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
        }
 
        // Dummy index renaming

-       r = rename_dummy_indices(r, dummy_indices, local_dummy_indices);
+       r = rename_dummy_indices<idx>(r, dummy_indices, local_dummy_indices);
+       r = rename_dummy_indices<varidx>(r, dummy_indices, local_dummy_indices);
+       r = rename_dummy_indices<spinidx>(r, dummy_indices, local_dummy_indices);

 
        // Product of indexed object with a scalar?
        if (is_exactly_a<mul>(r) && r.nops() == 2
 
        // Product of indexed object with a scalar?
        if (is_exactly_a<mul>(r) && r.nops() == 2


@@ -943,6 +984,17 @@ public:
        }
 };
 
        }
 };
 

+bool hasindex(const ex &x, const ex &sym)
+{      
+       if(is_a<idx>(x) && x.op(0)==sym)
+               return true;
+       else
+               for(size_t i=0; i<x.nops(); ++i)
+                       if(hasindex(x.op(i), sym))
+                               return true;
+       return false;
+}
+

 /** Simplify indexed expression, return list of free indices. */
 ex simplify_indexed(const ex & e, exvector & free_indices, exvector & dummy_indices, const scalar_products & sp)
 {
 /** Simplify indexed expression, return list of free indices. */
 ex simplify_indexed(const ex & e, exvector & free_indices, exvector & dummy_indices, const scalar_products & sp)
 {


@@ -971,7 +1023,10 @@ ex simplify_indexed(const ex & e, exvector & free_indices, exvector & dummy_indi
                }
 
                // Rename the dummy indices
                }
 
                // Rename the dummy indices

-               return rename_dummy_indices(e_expanded, dummy_indices, local_dummy_indices);
+               e_expanded = rename_dummy_indices<idx>(e_expanded, dummy_indices, local_dummy_indices);
+               e_expanded = rename_dummy_indices<varidx>(e_expanded, dummy_indices, local_dummy_indices);
+               e_expanded = rename_dummy_indices<spinidx>(e_expanded, dummy_indices, local_dummy_indices);
+               return e_expanded;

        }
 
        // Simplification of sum = sum of simplifications, check consistency of
        }
 
        // Simplification of sum = sum of simplifications, check consistency of


@@ -1015,18 +1070,19 @@ ex simplify_indexed(const ex & e, exvector & free_indices, exvector & dummy_indi
                if (num_terms_orig < 2 || dummy_indices.size() < 2)
                        return sum;
 
                if (num_terms_orig < 2 || dummy_indices.size() < 2)
                        return sum;
 

-               // Yes, construct vector of all dummy index symbols
-               exvector dummy_syms;
-               dummy_syms.reserve(dummy_indices.size());
-               for (exvector::const_iterator it = dummy_indices.begin(); it != dummy_indices.end(); ++it)
-                       dummy_syms.push_back(it->op(0));
-

                // Chop the sum into terms and symmetrize each one over the dummy
                // indices
                std::vector<terminfo> terms;
                for (size_t i=0; i<sum.nops(); i++) {
                        const ex & term = sum.op(i);
                // Chop the sum into terms and symmetrize each one over the dummy
                // indices
                std::vector<terminfo> terms;
                for (size_t i=0; i<sum.nops(); i++) {
                        const ex & term = sum.op(i);

-                       ex term_symm = symmetrize(term, dummy_syms);
+                       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);
+                       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);

                        if (term_symm.is_zero())
                                continue;
                        terms.push_back(terminfo(term, term_symm));
                        if (term_symm.is_zero())
                                continue;
                        terms.push_back(terminfo(term, term_symm));


@@ -1318,9 +1374,9 @@ exvector get_all_dummy_indices(const ex & e)
        return v;
 }
 
        return v;
 }
 

-ex rename_dummy_indices_uniquely(const ex & a, const ex & b)
+ex rename_dummy_indices_uniquely(const exvector & va, const exvector & vb, const ex & b)

 {
 {

-       exvector va = get_all_dummy_indices(a), vb = get_all_dummy_indices(b), common_indices;
+       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 b;
        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 b;


@@ -1330,15 +1386,25 @@ ex rename_dummy_indices_uniquely(const ex & a, const ex & b)
                new_indices.reserve(2*common_indices.size());
                exvector::const_iterator ip = common_indices.begin(), ipend = common_indices.end();
                while (ip != ipend) {
                new_indices.reserve(2*common_indices.size());
                exvector::const_iterator ip = common_indices.begin(), ipend = common_indices.end();
                while (ip != ipend) {

-                       if (is_a<varidx>(*ip)) {
-                               varidx mu((new symbol)->setflag(status_flags::dynallocated), ex_to<varidx>(*ip).get_dim(), ex_to<varidx>(*ip).is_covariant());
-                               old_indices.push_back(*ip);
-                               new_indices.push_back(mu);
+                       ex newsym=(new symbol)->setflag(status_flags::dynallocated);
+                       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);
+                       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);
+                       else
+                               newidx = (new idx(newsym, ex_to<idx>(*ip).get_dim()))
+                                       -> setflag(status_flags::dynallocated);
+                       old_indices.push_back(*ip);
+                       new_indices.push_back(newidx);
+                       if(is_a<varidx>(*ip)) {

                                old_indices.push_back(ex_to<varidx>(*ip).toggle_variance());
                                old_indices.push_back(ex_to<varidx>(*ip).toggle_variance());

-                               new_indices.push_back(mu.toggle_variance());
-                       } else {
-                               old_indices.push_back(*ip);
-                               new_indices.push_back(idx((new symbol)->setflag(status_flags::dynallocated), ex_to<varidx>(*ip).get_dim()));
+                               new_indices.push_back(ex_to<varidx>(newidx).toggle_variance());

                        }
                        ++ip;
                }
                        }
                        ++ip;
                }


@@ -1346,6 +1412,15 @@ ex rename_dummy_indices_uniquely(const ex & a, const ex & b)
        }
 }
 
        }
 }
 

+ex rename_dummy_indices_uniquely(const ex & a, const ex & b)
+{
+       exvector va = get_all_dummy_indices(a);
+       exvector vb = get_all_dummy_indices(b);
+       sort(va.begin(), va.end(), ex_is_less());
+       sort(vb.begin(), vb.end(), ex_is_less());
+       return rename_dummy_indices_uniquely(va, vb, b);
+}
+

 ex expand_dummy_sum(const ex & e, bool subs_idx)
 {
        ex e_expanded = e.expand();
 ex expand_dummy_sum(const ex & e, bool subs_idx)
 {
        ex e_expanded = e.expand();

diff --git a/ginac/indexed.h b/ginac/indexed.h
index ad32c2372fee25e55e0fef167de8be32000eeea0..1ce81945b432cdc40f0537a7699eb35441043eb8 100644 (file)
--- a/ginac/indexed.h
+++ b/ginac/indexed.h
@@ -153,7 +153,8 @@ protected:
        ex derivative(const symbol & s) const;
        ex thiscontainer(const exvector & v) const;
        ex thiscontainer(std::auto_ptr<exvector> vp) const;
        ex derivative(const symbol & s) const;
        ex thiscontainer(const exvector & v) const;
        ex thiscontainer(std::auto_ptr<exvector> vp) const;

-       unsigned return_type() const { return return_types::commutative; }
+       unsigned return_type() const;
+       unsigned return_type_tinfo() const { return op(0).return_type_tinfo(); }

        ex expand(unsigned options = 0) const;
 
        // new virtual functions which can be overridden by derived classes
        ex expand(unsigned options = 0) const;
 
        // new virtual functions which can be overridden by derived classes


@@ -256,6 +257,9 @@ exvector get_all_dummy_indices(const ex & e);
 /** Returns b with all dummy indices, which are common with a, renamed */
 ex rename_dummy_indices_uniquely(const ex & a, const ex & b);
 
 /** Returns b with all dummy indices, which are common with a, renamed */
 ex rename_dummy_indices_uniquely(const ex & a, const ex & b);
 

+/** Same as above, where va and vb contain the indices of a and b and are sorted */
+ex rename_dummy_indices_uniquely(const exvector & va, const exvector & vb, const ex & b);
+

 /** This function returns the given expression with expanded sums
  *  for all dummy index summations, where the dimensionality of 
  *  the dummy index is numeric.
 /** This function returns the given expression with expanded sums
  *  for all dummy index summations, where the dimensionality of 
  *  the dummy index is numeric.

diff --git a/ginac/symmetry.cpp b/ginac/symmetry.cpp
index f29c6a38d1445096baf62bb104bd1cf2aca35d0b..39b194a6933cd274cb45d7dfe8033186b478bbef 100644 (file)
--- a/ginac/symmetry.cpp
+++ b/ginac/symmetry.cpp
@@ -427,7 +427,7 @@ static ex symm(const ex & e, exvector::const_iterator first, exvector::const_ite
                lst new_lst;
                for (unsigned i=0; i<num; i++)
                        new_lst.append(orig_lst.op(iv[i]));
                lst new_lst;
                for (unsigned i=0; i<num; i++)
                        new_lst.append(orig_lst.op(iv[i]));

-               ex term = e.subs(orig_lst, new_lst, subs_options::no_pattern);
+               ex term = e.subs(orig_lst, new_lst, subs_options::no_pattern|subs_options::no_index_renaming);

                if (asymmetric) {
                        memcpy(iv2, iv, num * sizeof(unsigned));
                        term *= permutation_sign(iv2, iv2 + num);
                if (asymmetric) {
                        memcpy(iv2, iv, num * sizeof(unsigned));
                        term *= permutation_sign(iv2, iv2 + num);


@@ -468,7 +468,7 @@ ex symmetrize_cyclic(const ex & e, exvector::const_iterator first, exvector::con
        for (unsigned i=0; i<num-1; i++) {
                ex perm = new_lst.op(0);
                new_lst.remove_first().append(perm);
        for (unsigned i=0; i<num-1; i++) {
                ex perm = new_lst.op(0);
                new_lst.remove_first().append(perm);

-               sum += e.subs(orig_lst, new_lst, subs_options::no_pattern);
+               sum += e.subs(orig_lst, new_lst, subs_options::no_pattern|subs_options::no_index_renaming);

        }
        return sum / num;
 }
        }
        return sum / num;
 }

diff --git a/ginac/tensor.h b/ginac/tensor.h
index 8078b0937a130aec0a622a8ea1a85efb47587edd..0ebfb0f56d76b93633c8a635d30c20e2cf763982 100644 (file)
--- a/ginac/tensor.h
+++ b/ginac/tensor.h
@@ -65,6 +65,7 @@ public:
 
        // non-virtual functions in this class
 protected:
 
        // non-virtual functions in this class
 protected:

+       unsigned return_type() const { return return_types::commutative; }

        void do_print(const print_context & c, unsigned level) const;
        void do_print_latex(const print_latex & c, unsigned level) const;
 };
        void do_print(const print_context & c, unsigned level) const;
        void do_print_latex(const print_latex & c, unsigned level) const;
 };


@@ -84,6 +85,7 @@ public:
 
        // non-virtual functions in this class
 protected:
 
        // non-virtual functions in this class
 protected:

+       unsigned return_type() const { return return_types::commutative; }

        void do_print(const print_context & c, unsigned level) const;
 };
 
        void do_print(const print_context & c, unsigned level) const;
 };
 


@@ -106,6 +108,7 @@ public:
 
        // non-virtual functions in this class
 protected:
 
        // non-virtual functions in this class
 protected:

+       unsigned return_type() const { return return_types::commutative; }

        void do_print(const print_context & c, unsigned level) const;
        void do_print_latex(const print_latex & c, unsigned level) const;
 
        void do_print(const print_context & c, unsigned level) const;
        void do_print_latex(const print_latex & c, unsigned level) const;
 


@@ -153,6 +156,7 @@ public:
 
        // non-virtual functions in this class
 protected:
 
        // non-virtual functions in this class
 protected:

+       unsigned return_type() const { return return_types::commutative; }

        void do_print(const print_context & c, unsigned level) const;
        void do_print_latex(const print_latex & c, unsigned level) const;
 
        void do_print(const print_context & c, unsigned level) const;
        void do_print_latex(const print_latex & c, unsigned level) const;