- added symmetrize() and antisymmetrize() functions

[ginac.git] / ginac / indexed.cpp

diff --git a/ginac/indexed.cpp b/ginac/indexed.cpp
index f8d090307a526d3f69336f8a392deb63ba07e21a..20e54752a93d89c90e17607c953f0e62741fbacd 100644 (file)
--- a/ginac/indexed.cpp
+++ b/ginac/indexed.cpp
@@ -30,6 +30,7 @@
 #include "ncmul.h"
 #include "power.h"
 #include "lst.h"
+#include "inifcns.h"
 #include "print.h"
 #include "archive.h"
 #include "utils.h"
@@ -479,6 +480,17 @@ exvector indexed::get_dummy_indices(const indexed & other) const
        return dummy_indices;
 }
 
+bool indexed::has_dummy_index_for(const ex & i) const
+{
+       exvector::const_iterator it = seq.begin() + 1, itend = seq.end();
+       while (it != itend) {
+               if (is_dummy_pair(*it, i))
+                       return true;
+               it++;
+       }
+       return false;
+}
+
 exvector indexed::get_free_indices(void) const
 {
        exvector free_indices, dummy_indices;
@@ -632,6 +644,8 @@ try_again:
                if (!is_ex_of_type(*it1, indexed))
                        continue;
 
+               bool first_noncommutative = (it1->return_type() != return_types::commutative);
+
                // Indexed factor found, get free indices and look for contraction
                // candidates
                exvector free1, dummy1;
@@ -643,6 +657,8 @@ try_again:
                        if (!is_ex_of_type(*it2, indexed))
                                continue;
 
+                       bool second_noncommutative = (it2->return_type() != return_types::commutative);
+
                        // Find free indices of second factor and merge them with free
                        // indices of first factor
                        exvector un;
@@ -685,7 +701,7 @@ try_again:
                        }
                        if (contracted) {
 contraction_done:
-                               if (non_commutative
+                               if (first_noncommutative || second_noncommutative
                                 || 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)
                                 || is_ex_exactly_of_type(*it1, ncmul) || is_ex_exactly_of_type(*it2, ncmul)) {
@@ -695,7 +711,7 @@ contraction_done:
                                        // Non-commutative products are always re-expanded to give
                                        // simplify_ncmul() the chance to re-order and canonicalize
                                        // the product
-                                       ex r = (non_commutative ? ex(ncmul(v)) : ex(mul(v)));
+                                       ex r = (non_commutative ? ex(ncmul(v, true)) : ex(mul(v)));
                                        return simplify_indexed(r, free_indices, dummy_indices, sp);
                                }
 
@@ -729,7 +745,7 @@ contraction_done:
 
        ex r;
        if (something_changed)
-               r = non_commutative ? ex(ncmul(v)) : ex(mul(v));
+               r = non_commutative ? ex(ncmul(v, true)) : ex(mul(v));
        else
                r = e;
 
@@ -812,6 +828,16 @@ ex simplify_indexed(const ex & e, const scalar_products & sp)
        return simplify_indexed(e, free_indices, dummy_indices, sp);
 }
 
+ex symmetrize(const ex & e)
+{
+       return symmetrize(e, e.get_free_indices());
+}
+
+ex antisymmetrize(const ex & e)
+{
+       return antisymmetrize(e, e.get_free_indices());
+}
+
 //////////
 // helper classes
 //////////