- dummy index renaming works better
authorChristian Bauer <Christian.Bauer@uni-mainz.de>
Mon, 23 Jul 2001 21:08:46 +0000 (21:08 +0000)
committerChristian Bauer <Christian.Bauer@uni-mainz.de>
Mon, 23 Jul 2001 21:08:46 +0000 (21:08 +0000)
- implemented contractions of epsilon tensors

ginac/indexed.cpp
ginac/tensor.cpp
ginac/tensor.h

index c6a84e3..56fd0b6 100644 (file)
@@ -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
index 65c8b59..28a7d7b 100644 (file)
@@ -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
 //////////
index d868355..e30a91a 100644 (file)
@@ -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: