simplify_indexed() will re-expand the expression and re-run the

[ginac.git] / ginac / indexed.cpp

diff --git a/ginac/indexed.cpp b/ginac/indexed.cpp
index 0ab989e801baccc8a94d60a7df7306da1822a78f..833d46c94111db070e6acf0ee34c6e8120e5f9e1 100644 (file)
--- a/ginac/indexed.cpp
+++ b/ginac/indexed.cpp
@@ -350,6 +350,24 @@ ex indexed::eval(int level) const
        return base.bp->eval_indexed(*this);
 }
 
+int indexed::degree(const ex & s) const
+{
+       return is_equal(*s.bp) ? 1 : 0;
+}
+
+int indexed::ldegree(const ex & s) const
+{
+       return is_equal(*s.bp) ? 1 : 0;
+}
+
+ex indexed::coeff(const ex & s, int n) const
+{
+       if (is_equal(*s.bp))
+               return n==1 ? _ex1() : _ex0();
+       else
+               return n==0 ? ex(*this) : _ex0();
+}
+
 ex indexed::thisexprseq(const exvector & v) const
 {
        return indexed(symmetry, v);
@@ -587,28 +605,36 @@ try_again:
                if (!is_ex_of_type(*it1, indexed))
                        continue;
 
-               // Indexed factor found, look for contraction candidates
+               // Indexed factor found, get free indices and look for contraction
+               // candidates
+               exvector free1, dummy1;
+               find_free_and_dummy(ex_to_indexed(*it1).seq.begin() + 1, ex_to_indexed(*it1).seq.end(), free1, dummy1);
+
                exvector::iterator it2;
                for (it2 = it1 + 1; it2 != itend; it2++) {
 
                        if (!is_ex_of_type(*it2, indexed))
                                continue;
 
+                       // Find free indices of second factor and merge them with free
+                       // indices of first factor
+                       exvector un;
+                       find_free_and_dummy(ex_to_indexed(*it2).seq.begin() + 1, ex_to_indexed(*it2).seq.end(), un, dummy1);
+                       un.insert(un.end(), free1.begin(), free1.end());
+
                        // Check whether the two factors share dummy indices
-                       exvector un(ex_to_indexed(*it1).seq.begin() + 1, ex_to_indexed(*it1).seq.end());
-                       un.insert(un.end(), ex_to_indexed(*it2).seq.begin() + 1, ex_to_indexed(*it2).seq.end());
                        exvector free, dummy;
                        find_free_and_dummy(un, free, dummy);
                        if (dummy.size() == 0)
                                continue;
 
                        // At least one dummy index, is it a defined scalar product?
+                       bool contracted = false;
                        if (free.size() == 0) {
                                if (sp.is_defined(*it1, *it2)) {
                                        *it1 = sp.evaluate(*it1, *it2);
                                        *it2 = _ex1();
-                                       something_changed = true;
-                                       goto try_again;
+                                       goto contraction_done;
                                }
                        }
 
@@ -623,7 +649,7 @@ try_again:
                        }
 
                        // Try to contract the first one with the second one
-                       bool contracted = it1->op(0).bp->contract_with(it1, it2, v);
+                       contracted = it1->op(0).bp->contract_with(it1, it2, v);
                        if (!contracted) {
 
                                // That didn't work; maybe the second object knows how to
@@ -631,11 +657,20 @@ try_again:
                                contracted = it2->op(0).bp->contract_with(it2, it1, v);
                        }
                        if (contracted) {
-                               something_changed = true;
+contraction_done:
+                               if (is_ex_exactly_of_type(*it1, add) || is_ex_exactly_of_type(*it2, add)
+                                || is_ex_exactly_of_type(*it1, mul) || is_ex_exactly_of_type(*it2, mul)) {
+
+                                       // One of the factors became a sum or product:
+                                       // re-expand expression and run again
+                                       ex r = non_commutative ? ex(ncmul(v)) : ex(mul(v));
+                                       return simplify_indexed(r, free_indices, sp);
+                               }
 
                                // Both objects may have new indices now or they might
                                // even not be indexed objects any more, so we have to
                                // start over
+                               something_changed = true;
                                goto try_again;
                        }
                }
@@ -645,21 +680,16 @@ try_again:
        exvector un, dummy_indices;
        it1 = v.begin(); itend = v.end();
        while (it1 != itend) {
-               if (is_ex_of_type(*it1, indexed)) {
-                       const indexed & o = ex_to_indexed(*it1);
-                       un.insert(un.end(), o.seq.begin() + 1, o.seq.end());
-               }
+               exvector free_indices_of_factor = it1->get_free_indices();
+               un.insert(un.end(), free_indices_of_factor.begin(), free_indices_of_factor.end());
                it1++;
        }
        find_free_and_dummy(un, free_indices, dummy_indices);
 
        ex r;
-       if (something_changed) {
-               if (non_commutative)
-                       r = ncmul(v);
-               else
-                       r = mul(v);
-       } else
+       if (something_changed)
+               r = non_commutative ? ex(ncmul(v)) : ex(mul(v));
+       else
                r = e;
 
        // Product of indexed object with a scalar?