From: Christian Bauer <Christian.Bauer@uni-mainz.de>
Date: Sun, 27 May 2001 19:53:13 +0000 (+0000)
Subject: - added symmetrize() and antisymmetrize() functions
X-Git-Tag: release_0-9-0~32
X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=commitdiff_plain;h=0117bd6ef4af029934703940d59e1c70866937b0

- added symmetrize() and antisymmetrize() functions
- generalized permutation_sign() template
- moved shaker_sort() template to utils.h
---

diff --git a/NEWS b/NEWS
index d08e730e..ec456f49 100644
--- a/NEWS
+++ b/NEWS
@@ -10,6 +10,7 @@ This file records noteworthy changes.
   optionally containing wildcard objects. These are constructed with the
   call "wild(<unsigned>)" and are denoted as "$0", "$1" etc. in the output
   and in ginsh.
+* Added symmetrize() and antisymmetrize() functions.
 * Fixed possible crash when calling subs() on expressions with non-commutative
   products.
 
diff --git a/check/exam_indexed.cpp b/check/exam_indexed.cpp
index 10eb3d01..75a99e51 100644
--- a/check/exam_indexed.cpp
+++ b/check/exam_indexed.cpp
@@ -161,6 +161,15 @@ static unsigned symmetry_check(void)
 	    indexed(B, indexed::antisymmetric, k, l); // GiNaC 0.8.0 had a bug here
 	result += check_equal_simplify(e, e);
 
+	e = indexed(A, i, j);
+	result += check_equal(symmetrize(e) + antisymmetrize(e), e);
+	e = indexed(A, indexed::symmetric, i, j, k, l);
+	result += check_equal(symmetrize(e), e);
+	result += check_equal(antisymmetrize(e), 0);
+	e = indexed(A, indexed::antisymmetric, i, j, k, l);
+	result += check_equal(symmetrize(e), 0);
+	result += check_equal(antisymmetrize(e), e);
+
 	return result;
 }
 
diff --git a/ginac/clifford.cpp b/ginac/clifford.cpp
index e9b50e9a..b5dff072 100644
--- a/ginac/clifford.cpp
+++ b/ginac/clifford.cpp
@@ -446,7 +446,7 @@ ex dirac_trace(const ex & e, unsigned char rl, const ex & trONE)
 								iv.push_back(n);
 								v.push_back(e.op(n));
 							}
-							int sign = permutation_sign(iv);
+							int sign = permutation_sign(iv.begin(), iv.end());
 							result += sign * eps0123(idx1, idx2, idx3, idx4)
 							        * dirac_trace(ncmul(v, true), rl, trONE);
 						}
diff --git a/ginac/idx.cpp b/ginac/idx.cpp
index 7088c5bd..c72e478e 100644
--- a/ginac/idx.cpp
+++ b/ginac/idx.cpp
@@ -472,36 +472,6 @@ bool is_dummy_pair(const ex & e1, const ex & e2)
 	return is_dummy_pair(ex_to_idx(e1), ex_to_idx(e2));
 }
 
-// Shaker sort is sufficient for the expected small number of indices
-template <class It, class Cmp>
-inline void shaker_sort(It first, It last, Cmp comp)
-{
-	if (first == last)
-		return;
-	--last;
-	if (first == last)
-		return;
-	It flag = first;
-	do {
-		It i;
-		for (i=last; i>first; --i) {
-			if (comp(*i, i[-1])) {
-				iter_swap(i-1, i);
-				flag = i - 1;
-			}
-		}
-		++flag;
-		first = flag;
-		for (i=first; i<last; ++i) {
-			if (comp(i[1], *i)) {
-				iter_swap(i, i+1);
-				flag = i + 1;
-			}
-		}
-		last = flag - 1;
-	} while (first <= last);
-}
-
 void find_free_and_dummy(exvector::const_iterator it, exvector::const_iterator itend, exvector & out_free, exvector & out_dummy)
 {
 	out_free.clear();
diff --git a/ginac/indexed.cpp b/ginac/indexed.cpp
index 9ddff040..20e54752 100644
--- 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"
@@ -827,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
 //////////
diff --git a/ginac/indexed.h b/ginac/indexed.h
index 9f6af5fd..fd85b8db 100644
--- a/ginac/indexed.h
+++ b/ginac/indexed.h
@@ -263,6 +263,11 @@ ex simplify_indexed(const ex & e);
  *  @return simplified expression */
 ex simplify_indexed(const ex & e, const scalar_products & sp);
 
+/** Symmetrize expression over its free indices. */
+ex symmetrize(const ex & e);
+
+/** Antisymmetrize expression over its free indices. */
+ex antisymmetrize(const ex & e);
 
 } // namespace GiNaC
 
diff --git a/ginac/inifcns.cpp b/ginac/inifcns.cpp
index b55261eb..26672154 100644
--- a/ginac/inifcns.cpp
+++ b/ginac/inifcns.cpp
@@ -523,6 +523,71 @@ ex ncpow(const ex & basis, unsigned exponent)
 	return ncmul(v, true);
 }
 
+// Symmetrize/antisymmetrize over a vector of objects
+static ex symm(const ex & e, exvector::const_iterator first, exvector::const_iterator last, bool asymmetric)
+{
+	// Need at least 2 objects for this operation
+	int num = last - first;
+	if (num < 2)
+		return e;
+
+	// Sort object vector, transform it into a list, and make a copy so we
+	// will know which objects get substituted for which
+	exvector iv(first, last);
+	sort(iv.begin(), iv.end(), ex_is_less());
+	exlist iv_lst;
+	iv_lst.insert(iv_lst.begin(), iv.begin(), iv.end());
+	lst orig_lst(iv_lst);
+
+	// With n objects there are n! possible permutations
+	int num_perms = factorial(numeric(num)).to_int();
+
+	// Loop over all permutations (the first permutation, which is the
+	// identity, is unrolled)
+	ex sum = e;
+	int i = 1;
+	do {
+		next_permutation(iv_lst.begin(), iv_lst.end(), ex_is_less());
+		ex term = e.subs(orig_lst, lst(iv_lst));
+		if (asymmetric) {
+			exlist test_lst = iv_lst;
+			term *= permutation_sign(test_lst.begin(), test_lst.end(), ex_is_less());
+		}
+		sum += term;
+		i++;
+	} while (i < num_perms);
+
+	return sum / num_perms;
+}
+
+ex symmetrize(const ex & e, exvector::const_iterator first, exvector::const_iterator last)
+{
+	return symm(e, first, last, false);
+}
+
+ex antisymmetrize(const ex & e, exvector::const_iterator first, exvector::const_iterator last)
+{
+	return symm(e, first, last, true);
+}
+
+ex symmetrize(const ex & e, const lst & l)
+{
+	exvector v;
+	v.reserve(l.nops());
+	for (unsigned i=0; i<l.nops(); i++)
+		v.push_back(l.op(i));
+	return symm(e, v.begin(), v.end(), false);
+}
+
+ex antisymmetrize(const ex & e, const lst & l)
+{
+	exvector v;
+	v.reserve(l.nops());
+	for (unsigned i=0; i<l.nops(); i++)
+		v.push_back(l.op(i));
+	return symm(e, v.begin(), v.end(), true);
+}
+
 /** Force inclusion of functions from initcns_gamma and inifcns_zeta
  *  for static lib (so ginsh will see them). */
 unsigned force_include_tgamma = function_index_tgamma;
diff --git a/ginac/inifcns.h b/ginac/inifcns.h
index 985b50ff..6a7efbe6 100644
--- a/ginac/inifcns.h
+++ b/ginac/inifcns.h
@@ -136,6 +136,30 @@ ex lsolve(const ex &eqns, const ex &symbols);
 /** Power of non-commutative basis. */
 ex ncpow(const ex & basis, unsigned exponent);
 
+/** Symmetrize expression over a set of objects (symbols, indices). */
+ex symmetrize(const ex & e, exvector::const_iterator first, exvector::const_iterator last);
+
+/** Symmetrize expression over a set of objects (symbols, indices). */
+inline ex symmetrize(const ex & e, const exvector & v)
+{
+	return symmetrize(e, v.begin(), v.end());
+}
+
+/** Symmetrize expression over a list of objects (symbols, indices). */
+ex symmetrize(const ex & e, const lst & l);
+
+/** Antisymmetrize expression over a set of objects (symbols, indices). */
+ex antisymmetrize(const ex & e, exvector::const_iterator first, exvector::const_iterator last);
+
+/** Antisymmetrize expression over a set of objects (symbols, indices). */
+inline ex antisymmetrize(const ex & e, const exvector & v)
+{
+	return antisymmetrize(e, v.begin(), v.end());
+}
+
+/** Antisymmetrize expression over a list of objects (symbols, indices). */
+ex antisymmetrize(const ex & e, const lst & l);
+
 inline bool is_order_function(const ex & e)
 {
 	return is_ex_the_function(e, Order);
diff --git a/ginac/matrix.cpp b/ginac/matrix.cpp
index 65b12549..99083c06 100644
--- a/ginac/matrix.cpp
+++ b/ginac/matrix.cpp
@@ -756,7 +756,7 @@ ex matrix::determinant(unsigned algo) const
 			std::vector<unsigned> pre_sort;
 			for (std::vector<uintpair>::iterator i=c_zeros.begin(); i!=c_zeros.end(); ++i)
 				pre_sort.push_back(i->second);
-			int sign = permutation_sign(pre_sort);
+			int sign = permutation_sign(pre_sort.begin(), pre_sort.end());
 			exvector result(row*col);  // represents sorted matrix
 			unsigned c = 0;
 			for (std::vector<unsigned>::iterator i=pre_sort.begin();
diff --git a/ginac/tensor.cpp b/ginac/tensor.cpp
index 23be5b0b..f4e49451 100644
--- a/ginac/tensor.cpp
+++ b/ginac/tensor.cpp
@@ -305,7 +305,7 @@ ex tensepsilon::eval_indexed(const basic & i) const
 		v.reserve(i.nops() - 1);
 		for (unsigned j=1; j<i.nops(); j++)
 			v.push_back(ex_to_numeric(ex_to_idx(i.op(j)).get_value()).to_int());
-		int sign = permutation_sign(v);
+		int sign = permutation_sign(v.begin(), v.end());
 
 		// In a Minkowski space, check for covariant indices
 		if (minkowski) {
diff --git a/ginac/utils.h b/ginac/utils.h
index 512c97c3..efdeaa2b 100644
--- a/ginac/utils.h
+++ b/ginac/utils.h
@@ -127,24 +127,106 @@ inline unsigned golden_ratio_hash(unsigned n)
 #endif
 }
 
-// Compute the sign of a permutation of a vector of things.
-template <typename T>
-int permutation_sign(std::vector<T> s)
+/* Compute the sign of a permutation of a container, with and without an
+   explicitly supplied comparison function. The containers gets modified
+   during the operation. */
+template <class It>
+int permutation_sign(It first, It last)
 {
-	if (s.size() < 2)
+	if (first == last)
 		return 0;
-	int sigma = 1;
-	for (typename std::vector<T>::iterator i=s.begin(); i!=s.end()-1; ++i) {
-		for (typename std::vector<T>::iterator j=i+1; j!=s.end(); ++j) {
-			if (*i == *j)
-				return 0;
-			if (*i > *j) {
-				iter_swap(i,j);
-				sigma = -sigma;
+	It i = first;
+	++i;
+	if (i == last)
+		return 0;
+	i = first;
+	It next_to_last = last;
+	--next_to_last;
+
+	int sign = 1;
+	while (i != next_to_last) {
+		It j = i;
+		++j;
+		while (j != last) {
+			if (!(*i < *j)) {
+				if (!(*j < *i))
+					return 0;
+				iter_swap(i, j);
+				sign = -sign;
+			}
+			++j;
+		}
+		++i;
+	}
+	return sign;
+}
+
+/** Compute the sign of a permutation of a container */
+template <class It, class Cmp>
+int permutation_sign(It first, It last, Cmp comp)
+{
+	if (first == last)
+		return 0;
+	It i = first;
+	++i;
+	if (i == last)
+		return 0;
+	i = first;
+	It next_to_last = last;
+	--next_to_last;
+
+	int sign = 1;
+	while (i != next_to_last) {
+		It j = i;
+		++j;
+		while (j != last) {
+			if (!comp(*i, *j)) {
+				if (!comp(*j, *i))
+					return 0;
+				iter_swap(i, j);
+				sign = -sign;
 			}
+			++j;
 		}
+		++i;
 	}
-	return sigma;
+	return sign;
+}
+
+/* Implementation of shaker sort, only compares adjacent elements. */
+template <class It, class Cmp>
+inline void shaker_sort(It first, It last, Cmp comp)
+{
+	if (first == last)
+		return;
+	--last;
+	if (first == last)
+		return;
+	It flag = first;
+	do {
+		It i = last, other = last;
+		--other;
+		while (i > first) {
+			if (comp(*i, *other)) {
+				iter_swap(other, i);
+				flag = other;
+			}
+			--i; --other;
+		}
+		++flag;
+		first = flag;
+		i = first; other = first;
+		++other;
+		while (i < last) {
+			if (comp(*other, *i)) {
+				iter_swap(i, other);
+				flag = other;
+			}
+			++i; ++other;
+		}
+		last = flag;
+		--last;
+	} while (first <= last);
 }
 
 /* Function objects for STL sort() etc. */