From c4c7b40ea71dd20ecc8942f9cab50c4bb9766624 Mon Sep 17 00:00:00 2001
From: Christian Bauer <Christian.Bauer@uni-mainz.de>
Date: Mon, 23 Jul 2001 21:08:46 +0000
Subject: [PATCH 1/1] - dummy index renaming works better - implemented
 contractions of epsilon tensors

---
 ginac/indexed.cpp | 44 +++++++++++++++++++++++++++-----------------
 ginac/tensor.cpp  | 26 ++++++++++++++++++++++++++
 ginac/tensor.h    |  1 +
 3 files changed, 54 insertions(+), 17 deletions(-)

diff --git a/ginac/indexed.cpp b/ginac/indexed.cpp
index c6a84e31..56fd0b6b 100644
--- 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
+		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) {
@@ -555,26 +556,35 @@ static ex rename_dummy_indices(const ex & e, exvector & global_dummy_indices, ex
 			}
 			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;
-	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
diff --git a/ginac/tensor.cpp b/ginac/tensor.cpp
index 65c8b595..28a7d7bf 100644
--- a/ginac/tensor.cpp
+++ b/ginac/tensor.cpp
@@ -30,6 +30,7 @@
 #include "relational.h"
 #include "lst.h"
 #include "numeric.h"
+#include "matrix.h"
 #include "print.h"
 #include "archive.h"
 #include "utils.h"
@@ -493,6 +494,31 @@ again:
 	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
 //////////
diff --git a/ginac/tensor.h b/ginac/tensor.h
index d868355b..e30a91a4 100644
--- 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;
+	bool contract_with(exvector::iterator self, exvector::iterator other, exvector & v) const;
 
 	// member variables
 private:
-- 
2.44.0