- dummy index renaming works better

- implemented contractions of epsilon tensors

diff --git a/ginac/indexed.cpp b/ginac/indexed.cpp
index c6a84e31ac70e6346091651bd0eb61f1a70fef5d..56fd0b6b4080ebaafa72845978bfea39d5198ae9 100644 (file)
--- a/ginac/indexed.cpp
+++ b/ginac/indexed.cpp
@@ -545,6 +545,7 @@ static ex rename_dummy_indices(const ex & e, exvector & global_dummy_indices, ex
 
                // More local indices than we encountered before, add the new ones
                // to the global set
 
                // More local indices than we encountered before, add the new ones
                // to the global set

+               int old_global_size = global_size;

                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) {


@@ -555,26 +556,35 @@ static ex rename_dummy_indices(const ex & e, exvector & global_dummy_indices, ex
                        }
                        it++;
                }
                        }
                        it++;
                }

-       }
+               shaker_sort(global_dummy_indices.begin(), global_dummy_indices.end(), ex_is_less(), ex_swap());

 
 

-       // Replace index symbols in expression
-       GINAC_ASSERT(local_size <= global_size);
-       bool all_equal = true;
-       lst local_syms, global_syms;
-       for (unsigned i=0; i<local_size; i++) {
-               ex loc_sym = local_dummy_indices[i].op(0);
-               ex glob_sym = global_dummy_indices[i].op(0);
-               if (!loc_sym.is_equal(glob_sym)
-                && ex_to<idx>(local_dummy_indices[i]).get_dim().is_equal(ex_to<idx>(global_dummy_indices[i]).get_dim())) {
-                       all_equal = false;
-                       local_syms.append(loc_sym);
-                       global_syms.append(glob_sym);
-               }
+               // If this is the first set of local indices, do nothing
+               if (old_global_size == 0)
+                       return e;

        }
        }

-       if (all_equal)
+       GINAC_ASSERT(local_size <= global_size);
+
+       // Construct lists of index symbols
+       exlist local_syms, global_syms;
+       for (unsigned i=0; i<local_size; i++)
+               local_syms.push_back(local_dummy_indices[i].op(0));
+       shaker_sort(local_syms.begin(), local_syms.end(), ex_is_less(), ex_swap());
+       for (unsigned i=0; i<global_size; i++)
+               global_syms.push_back(global_dummy_indices[i].op(0));
+
+       // Remove common indices
+       exlist local_uniq, global_uniq;
+       set_difference(local_syms.begin(), local_syms.end(), global_syms.begin(), global_syms.end(), std::back_insert_iterator<exlist>(local_uniq), ex_is_less());
+       set_difference(global_syms.begin(), global_syms.end(), local_syms.begin(), local_syms.end(), std::back_insert_iterator<exlist>(global_uniq), ex_is_less());
+
+       // Replace remaining non-common local index symbols by global ones
+       if (local_uniq.empty())

                return e;
                return e;

-       else
-               return e.subs(local_syms, global_syms);
+       else {
+               while (global_uniq.size() > local_uniq.size())
+                       global_uniq.pop_back();
+               return e.subs(lst(local_uniq), lst(global_uniq));
+       }

 }
 
 /** Simplify product of indexed expressions (commutative, noncommutative and
 }
 
 /** Simplify product of indexed expressions (commutative, noncommutative and

diff --git a/ginac/tensor.cpp b/ginac/tensor.cpp
index 65c8b59577c6e2cf0b20ebca8ad5e1fd9dee11dd..28a7d7bf61b7d8831ce8c3c312b4027668c49b73 100644 (file)
--- a/ginac/tensor.cpp
+++ b/ginac/tensor.cpp
@@ -30,6 +30,7 @@
 #include "relational.h"
 #include "lst.h"
 #include "numeric.h"
 #include "relational.h"
 #include "lst.h"
 #include "numeric.h"

+#include "matrix.h"

 #include "print.h"
 #include "archive.h"
 #include "utils.h"
 #include "print.h"
 #include "archive.h"
 #include "utils.h"


@@ -493,6 +494,31 @@ again:
        return false;
 }
 
        return false;
 }
 

+/** Contraction of epsilon tensor with something else. */
+bool tensepsilon::contract_with(exvector::iterator self, exvector::iterator other, exvector & v) const
+{
+       GINAC_ASSERT(is_ex_of_type(*self, indexed));
+       GINAC_ASSERT(is_ex_of_type(*other, indexed));
+       GINAC_ASSERT(is_ex_of_type(self->op(0), spinmetric));
+       unsigned num = self->nops() - 1;
+
+       if (is_ex_exactly_of_type(other->op(0), tensepsilon) && num+1 == other->nops()) {
+
+               // Contraction of two epsilon tensors is a determinant
+               ex dim = ex_to<idx>(self->op(1)).get_dim();
+               matrix M(num, num);
+               for (int i=0; i<num; i++)
+                       for (int j=0; j<num; j++)
+                               M(i, j) = delta_tensor(self->op(i+1), other->op(j+1));
+               int sign = minkowski ? -1 : 1;
+               *self = sign * M.determinant().simplify_indexed();
+               *other = _ex1();
+               return true;
+       }
+
+       return false;
+}
+

 //////////
 // global functions
 //////////
 //////////
 // global functions
 //////////

diff --git a/ginac/tensor.h b/ginac/tensor.h
index d868355b16b71ee07e9102c1b1995dd5e7330792..e30a91a4dcf41ec84070f7ceeb0fca6b5f9fe949 100644 (file)
--- a/ginac/tensor.h
+++ b/ginac/tensor.h
@@ -128,6 +128,7 @@ public:
 public:
        void print(const print_context & c, unsigned level = 0) const;
        ex eval_indexed(const basic & i) const;
 public:
        void print(const print_context & c, unsigned level = 0) const;
        ex eval_indexed(const basic & i) const;

+       bool contract_with(exvector::iterator self, exvector::iterator other, exvector & v) const;

 
        // member variables
 private:
 
        // member variables
 private: