Fixing bug with respect to canonicalizing the variance of dummy indices in

the presence of a cyclic symmetry of indices.

diff --git a/ginac/indexed.cpp b/ginac/indexed.cpp
index bd41f44b0b38a271d243564c945d674bea9acaa3..d772b3f4743bd85ff8c127110f888c34801c0c6c 100644 (file)
--- a/ginac/indexed.cpp
+++ b/ginac/indexed.cpp
@@ -639,6 +639,56 @@ bool reposition_dummy_indices(ex & e, exvector & variant_dummy_indices, exvector
 {
        bool something_changed = false;
 
 {
        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;
+                                               }
+                                       }
+                                       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;
+                       }
+               }
+               e = optimal_e;
+       }
+

        // 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;
        // 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;

diff --git a/ginac/symmetry.cpp b/ginac/symmetry.cpp
index 9cee0ffb37b24ec718cabf91c3897307ff3c9fab..e48e269851e0846442bbbc1975928348e8d24271 100644 (file)
--- a/ginac/symmetry.cpp
+++ b/ginac/symmetry.cpp
@@ -266,6 +266,18 @@ void symmetry::do_print_tree(const print_tree & c, unsigned level) const
 // non-virtual functions in this class
 //////////
 
 // non-virtual functions in this class
 //////////
 

+bool symmetry::has_cyclic() const
+{
+       if (type == cyclic)
+               return true;
+
+       for (exvector::const_iterator i=children.begin(); i!=children.end(); ++i)
+               if (ex_to<symmetry>(*i).has_cyclic())
+                       return true;
+
+       return false;
+}
+

 symmetry &symmetry::add(const symmetry &c)
 {
        // All children must have the same number of indices
 symmetry &symmetry::add(const symmetry &c)
 {
        // All children must have the same number of indices

diff --git a/ginac/symmetry.h b/ginac/symmetry.h
index 4168013ed4d391e438154dae129107d4f702c535..80a2386d791048bf7da8fec7b4ad13c64d343cd5 100644 (file)
--- a/ginac/symmetry.h
+++ b/ginac/symmetry.h
@@ -80,6 +80,8 @@ public:
 
        /** Check whether this node actually represents any kind of symmetry. */
        bool has_symmetry() const {return type != none || !children.empty(); }
 
        /** Check whether this node actually represents any kind of symmetry. */
        bool has_symmetry() const {return type != none || !children.empty(); }

+       /** Check whether this node involves a cyclic symmetry. */
+       bool has_cyclic() const;

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