- dummy index renaming works better

[ginac.git] / ginac / tensor.cpp

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