- added symmetrize() and antisymmetrize() functions
authorChristian Bauer <Christian.Bauer@uni-mainz.de>
Sun, 27 May 2001 19:53:13 +0000 (19:53 +0000)
committerChristian Bauer <Christian.Bauer@uni-mainz.de>
Sun, 27 May 2001 19:53:13 +0000 (19:53 +0000)
- generalized permutation_sign() template
- moved shaker_sort() template to utils.h

NEWS
check/exam_indexed.cpp
ginac/clifford.cpp
ginac/idx.cpp
ginac/indexed.cpp
ginac/indexed.h
ginac/inifcns.cpp
ginac/inifcns.h
ginac/matrix.cpp
ginac/tensor.cpp
ginac/utils.h

diff --git a/NEWS b/NEWS
index d08e730e48711ef795f7a6e337820f202b3dcf0b..ec456f49047f1cf655a3605ab0f7d74b4ba71ff3 100644 (file)
--- 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.
 
index 10eb3d01e0b65b5dba44a18f4013cd264c918490..75a99e513a509410e4fb7c6826bf5d86f451deb9 100644 (file)
@@ -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;
 }
 
index e9b50e9a13c179feef05187cceb70e8ed4de1061..b5dff072b6601d2bb2893d7149e29eac50cdcf94 100644 (file)
@@ -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);
                                                }
index 7088c5bd674ae6584c76b14ce2669448e38e6e2f..c72e478e305bc02398394aa7486ceafe3d1d9150 100644 (file)
@@ -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();
index 9ddff040fd509c3d4f7773df99fcd7db7d34d8bf..20e54752a93d89c90e17607c953f0e62741fbacd 100644 (file)
@@ -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
 //////////
index 9f6af5fd1dd1001846da0ae5c3c322b584f6802e..fd85b8db25cfa055805e78a2547d63b40229b273 100644 (file)
@@ -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
 
index b55261eb1a0072939b0a1bb43a2a317d681778ec..2667215416129be392b21891aa2d16ed6bf89139 100644 (file)
@@ -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;
index 985b50ffe33cc45bb21ff43b4bea8df4538dfd09..6a7efbe610a24930690b1462eb42f02427585f9a 100644 (file)
@@ -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);
index 65b1254936e0a378eb5ba1760ee9709a3a52c409..99083c060be15496d74c849266fa5ce18e1624ff 100644 (file)
@@ -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();
index 23be5b0bed201caf8c465dca33be7bcaa1ef746f..f4e49451d907a70cbb541973984035cf872c6b1e 100644 (file)
@@ -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) {
index 512c97c34d93457435e7975fd876ee285d8d8c35..efdeaa2b7fe5cee59435176dcb32d216286f27c9 100644 (file)
@@ -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. */