Fixed bug in expanding expressions containing dummy indices. [V.Kisil]

diff --git a/ginac/indexed.cpp b/ginac/indexed.cpp
index 6d74b09a855321cb80aa11057c86c0c6b657833a..ef904f7222bfc798420dee5cdb6bed6373426a26 100644 (file)
--- a/ginac/indexed.cpp
+++ b/ginac/indexed.cpp
@@ -35,8 +35,10 @@
 #include "operators.h"
 #include "lst.h"
 #include "archive.h"
+#include "symbol.h"
 #include "utils.h"
 #include "integral.h"
+#include "matrix.h"
 
 namespace GiNaC {
 
@@ -669,16 +671,15 @@ struct ex_base_is_less : public std::binary_function<ex, ex, bool> {
        }
 };
 
-/** 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)
+/* An auxiliary function used by simplify_indexed() and expand_dummy_sum() 
+ * It returns an exvector of factors from the supplied product */
+static void product_to_exvector(const ex & e, exvector & v, bool & non_commutative)
 {
        // Remember whether the product was commutative or noncommutative
        // (because we chop it into factors and need to reassemble later)
-       bool non_commutative = is_exactly_a<ncmul>(e);
+       non_commutative = is_exactly_a<ncmul>(e);
 
        // Collect factors in an exvector, store squares twice
-       exvector v;
        v.reserve(e.nops() * 2);
 
        if (is_exactly_a<power>(e)) {
@@ -691,7 +692,7 @@ ex simplify_indexed_product(const ex & e, exvector & free_indices, exvector & du
                        ex f = e.op(i);
                        if (is_exactly_a<power>(f) && f.op(1).is_equal(_ex2)) {
                                v.push_back(f.op(0));
-                   v.push_back(f.op(0));
+                               v.push_back(f.op(0));
                        } else if (is_exactly_a<ncmul>(f)) {
                                // Noncommutative factor found, split it as well
                                non_commutative = true; // everything becomes noncommutative, ncmul will sort out the commutative factors later
@@ -701,6 +702,19 @@ ex simplify_indexed_product(const ex & e, exvector & free_indices, exvector & du
                                v.push_back(f);
                }
        }
+}
+
+/** 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)
+{
+       // Collect factors in an exvector
+       exvector v;
+
+       // Remember whether the product was commutative or noncommutative
+       // (because we chop it into factors and need to reassemble later)
+       bool non_commutative;
+       product_to_exvector(e, v, non_commutative);
 
        // Perform contractions
        bool something_changed = false;
@@ -1278,4 +1292,107 @@ void scalar_products::debugprint() const
        }
 }
 
+/** Returns all dummy indices from the exvector */
+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();
+       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;
+                       while (ip1 != ipend) {
+                               if (is_a<indexed>(*ip1)) {
+                                       v1 = ex_to<indexed>(*ip).get_dummy_indices(ex_to<indexed>(*ip1));
+                                       v.insert(v.end(), v1.begin(), v1.end());
+                               }
+                               ++ip1;
+                       }
+               }
+               ++ip;
+       }
+       return v;
+}
+
+ex rename_dummy_indices_uniquely(const ex & a, const ex & b)
+{
+       exvector va = get_all_dummy_indices(a), vb = get_all_dummy_indices(b), 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;
+       } 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) {
+                       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);
+                               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()));
+                       }
+                       ++ip;
+               }
+               return b.subs(lst(old_indices.begin(), old_indices.end()), lst(new_indices.begin(), new_indices.end()), subs_options::no_pattern);
+       }
+}
+
+ex expand_dummy_sum(const ex & e, bool subs_idx)
+{
+       ex e_expanded = e.expand();
+       pointer_to_map_function_1arg<bool> fcn(expand_dummy_sum, subs_idx);
+       if (is_a<add>(e_expanded) || is_a<lst>(e_expanded) || is_a<matrix>(e_expanded)) {
+               return e_expanded.map(fcn);
+       } else if (is_a<ncmul>(e_expanded) || is_a<mul>(e_expanded) || is_a<power>(e_expanded)) {
+               exvector v = get_all_dummy_indices(e_expanded);
+               exvector::const_iterator it = v.begin(), itend = v.end();
+               while (it != itend) {
+                       varidx nu = ex_to<varidx>(*it);
+                       if (nu.is_dim_numeric()) {
+                               ex en = 0;
+                               for (int i=0; i < ex_to<numeric>(nu.get_dim()).to_int(); i++) {
+                                       if (is_a<varidx>(nu) && !subs_idx) {
+                                               en += e_expanded.subs(lst(nu == varidx(i, nu.get_dim(), true), nu.toggle_variance() == varidx(i, nu.get_dim())));
+                                       } else {
+                                               en += e_expanded.subs(lst(nu == idx(i, nu.get_dim()), nu.toggle_variance() == idx(i, nu.get_dim())));
+                                       }
+                               }
+                               return expand_dummy_sum(en, subs_idx);
+                       } 
+                       ++it;
+               }
+               return e;
+       } else if (is_a<indexed>(e_expanded)) {
+               exvector v = ex_to<indexed>(e_expanded).get_dummy_indices();
+               exvector::const_iterator it = v.begin(), itend = v.end();
+               while (it != itend) {
+                       varidx nu = ex_to<varidx>(*it);
+                       if (nu.is_dim_numeric()) {
+                               ex en = 0;
+                               for (int i=0; i < ex_to<numeric>(nu.get_dim()).to_int(); i++) {
+                                       if (is_a<varidx>(nu) && !subs_idx) {
+                                               en += e_expanded.subs(lst(nu == varidx(i, nu.get_dim(), true), nu.toggle_variance() == varidx(i, nu.get_dim())));
+                                       } else {
+                                               en += e_expanded.subs(lst(nu == idx(i, nu.get_dim()), nu.toggle_variance() == idx(i, nu.get_dim())));
+                                       }
+                               }
+                               return expand_dummy_sum(en, subs_idx);
+                       } 
+                       ++it;
+               }
+               return e;
+       } else {
+               return e;
+       }
+}
+
 } // namespace GiNaC

diff --git a/ginac/indexed.h b/ginac/indexed.h
index 7329e87db1c7b64bda0eb3205cbaa75126b5e0b2..ad32c2372fee25e55e0fef167de8be32000eeea0 100644 (file)
--- a/ginac/indexed.h
+++ b/ginac/indexed.h
@@ -250,6 +250,23 @@ template<> inline bool is_exactly_a<indexed>(const basic & obj)
        return obj.tinfo()==TINFO_indexed;
 }
 
+/** Returns all dummy indices from the expression */
+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);
+
+/** This function returns the given expression with expanded sums
+ *  for all dummy index summations, where the dimensionality of 
+ *  the dummy index is numeric.
+ *  Optionally all indices with a variance will be substituted by 
+ *  indices with the corresponding numeric values without variance.
+ *
+ *  @param e the given expression
+ *  @param subs_idx indicates if variance of dummy indixes should be neglected
+ */
+ex expand_dummy_sum(const ex & e, bool subs_idx = false);
+
 } // namespace GiNaC
 
 #endif // ndef __GINAC_INDEXED_H__

diff --git a/ginac/mul.cpp b/ginac/mul.cpp
index 924493d272f1956f6d4d755f55f5b01203df5c48..fea5b346ce5ff3c3d9e922fb845166570cbd26b9 100644 (file)
--- a/ginac/mul.cpp
+++ b/ginac/mul.cpp
@@ -30,6 +30,7 @@
 #include "power.h"
 #include "operators.h"
 #include "matrix.h"
+#include "indexed.h"
 #include "lst.h"
 #include "archive.h"
 #include "utils.h"
@@ -904,11 +905,12 @@ ex mul::expand(unsigned options) const
                                        for (epvector::const_iterator i2=add2begin; i2!=add2end; ++i2) {
                                                // Don't push_back expairs which might have a rest that evaluates to a numeric,
                                                // since that would violate an invariant of expairseq:
-                                               const ex rest = (new mul(i1->rest, i2->rest))->setflag(status_flags::dynallocated);
-                                               if (is_exactly_a<numeric>(rest))
+                                               const ex rest = (new mul(i1->rest, rename_dummy_indices_uniquely(i1->rest, i2->rest)))->setflag(status_flags::dynallocated);
+                                               if (is_exactly_a<numeric>(rest)) {
                                                        oc += ex_to<numeric>(rest).mul(ex_to<numeric>(i1->coeff).mul(ex_to<numeric>(i2->coeff)));
-                                               else
+                                               } else {
                                                        distrseq.push_back(expair(rest, ex_to<numeric>(i1->coeff).mul_dyn(ex_to<numeric>(i2->coeff))));
+                                               }
                                        }
                                        tmp_accu += (new add(distrseq, oc))->setflag(status_flags::dynallocated);
                                }
@@ -934,11 +936,11 @@ ex mul::expand(unsigned options) const
 
                for (size_t i=0; i<n; ++i) {
                        epvector factors = non_adds;
-                       factors.push_back(split_ex_to_pair(last_expanded.op(i)));
+                       factors.push_back(split_ex_to_pair(rename_dummy_indices_uniquely(mul(non_adds), last_expanded.op(i))));
                        ex term = (new mul(factors, overall_coeff))->setflag(status_flags::dynallocated);
-                       if (can_be_further_expanded(term))
+                       if (can_be_further_expanded(term)) {
                                distrseq.push_back(term.expand());
-                       else {
+                       } else {
                                if (options == 0)
                                        ex_to<basic>(term).setflag(status_flags::expanded);
                                distrseq.push_back(term);

diff --git a/ginac/ncmul.cpp b/ginac/ncmul.cpp
index 4424253ce4369c9a1ebb5c271d12e500ca08638b..eeadf4e7feb4fa6acf7b48db84062b049e502526 100644 (file)
--- a/ginac/ncmul.cpp
+++ b/ginac/ncmul.cpp
@@ -30,6 +30,7 @@
 #include "mul.h"
 #include "matrix.h"
 #include "archive.h"
+#include "indexed.h"
 #include "utils.h"
 
 namespace GiNaC {
@@ -167,10 +168,25 @@ ex ncmul::expand(unsigned options) const
 
        intvector k(number_of_adds);
 
+       /* Rename indices in the static members of the product */
+       exvector expanded_seq_mod;
+       size_t j = 0;
+       size_t i;
+       for (i=0; i<expanded_seq.size(); i++) {
+               if (i == positions_of_adds[j]) {
+                       expanded_seq_mod.push_back(_ex1);
+                       j++;
+               } else {
+                       expanded_seq_mod.push_back(rename_dummy_indices_uniquely(ncmul(expanded_seq_mod), expanded_seq[i]));
+               }
+       }
+
        while (true) {
-               exvector term = expanded_seq;
-               for (size_t i=0; i<number_of_adds; i++)
-                       term[positions_of_adds[i]] = expanded_seq[positions_of_adds[i]].op(k[i]);
+               exvector term = expanded_seq_mod;
+               for (i=0; i<number_of_adds; i++) {
+                       term[positions_of_adds[i]] = rename_dummy_indices_uniquely(ncmul(term), expanded_seq[positions_of_adds[i]].op(k[i]));
+               }
+
                distrseq.push_back((new ncmul(term, true))->
                                    setflag(status_flags::dynallocated | (options == 0 ? status_flags::expanded : 0)));
 

diff --git a/ginac/power.cpp b/ginac/power.cpp
index bd276c56e454f035eb6043d1c9bf8b30ce769f72..7929768541fd30b9f90fd5693f9776793f9213dc 100644 (file)
--- a/ginac/power.cpp
+++ b/ginac/power.cpp
@@ -853,8 +853,17 @@ ex power::expand_mul(const mul & m, const numeric & n, unsigned options, bool fr
 {
        GINAC_ASSERT(n.is_integer());
 
-       if (n.is_zero())
+       if (n.is_zero()) {
                return _ex1;
+       }
+
+       // Leave it to multiplication since dummy indices have to be renamed
+       if (get_all_dummy_indices(m).size() > 0) {
+               ex result = m;
+               for (int i=1; i < n.to_int(); i++)
+                       result *= rename_dummy_indices_uniquely(m,m);
+               return result;
+       }
 
        epvector distrseq;
        distrseq.reserve(m.seq.size());