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.
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;
}
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);
}
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();
#include "ncmul.h"
#include "power.h"
#include "lst.h"
+#include "inifcns.h"
#include "print.h"
#include "archive.h"
#include "utils.h"
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
//////////
* @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
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;
/** 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);
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();
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) {
#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. */
git://www.ginac.de / ginac.git / commitdiff