- dirac_trace() is twice as fast

- permutation_sign() uses shaker sort
- shaker_sort() doesn't require less-than comparable iterators any more

diff --git a/ginac/clifford.cpp b/ginac/clifford.cpp
index b5dff072b6601d2bb2893d7149e29eac50cdcf94..3af4fcae1eab2b374e900d58192cf107a3406974 100644 (file)
--- a/ginac/clifford.cpp
+++ b/ginac/clifford.cpp
@@ -345,6 +345,42 @@ static bool is_clifford_tinfo(unsigned ti)
        return (ti & ~0xff) == TINFO_clifford;
 }
 
+/** Take trace of a string of an even number of Dirac gammas given a vector
+ *  of indices. */
+static ex trace_string(exvector::const_iterator ix, unsigned num)
+{
+       // Tr gamma.mu gamma.nu = 4 g.mu.nu
+       if (num == 2)
+               return lorentz_g(ix[0], ix[1]);
+
+       // Tr gamma.mu gamma.nu gamma.rho gamma.sig = 4 (g.mu.nu g.rho.sig + g.nu.rho g.mu.sig - g.mu.rho g.nu.sig
+       else if (num == 4)
+               return lorentz_g(ix[0], ix[1]) * lorentz_g(ix[2], ix[3])
+                    + lorentz_g(ix[1], ix[2]) * lorentz_g(ix[0], ix[3])
+                    - lorentz_g(ix[0], ix[2]) * lorentz_g(ix[1], ix[3]);
+
+       // Traces of 6 or more gammas are computed recursively:
+       // Tr gamma.mu1 gamma.mu2 ... gamma.mun =
+       //   + g.mu1.mu2 * Tr gamma.mu3 ... gamma.mun
+       //   - g.mu1.mu3 * Tr gamma.mu2 gamma.mu4 ... gamma.mun
+       //   + g.mu1.mu4 * Tr gamma.mu3 gamma.mu3 gamma.mu5 ... gamma.mun
+       //   - ...
+       //   + g.mu1.mun * Tr gamma.mu2 ... gamma.mu(n-1)
+       exvector v(num - 2);
+       int sign = 1;
+       ex result;
+       for (int i=1; i<num; i++) {
+               for (int n=1, j=0; n<num; n++) {
+                       if (n == i)
+                               continue;
+                       v[j++] = ix[n];
+               }
+               result += sign * lorentz_g(ix[0], ix[i]) * trace_string(v.begin(), num-2);
+               sign = -sign;
+       }
+       return result;
+}
+
 ex dirac_trace(const ex & e, unsigned char rl, const ex & trONE)
 {
        if (is_ex_of_type(e, clifford)) {
@@ -402,58 +438,41 @@ ex dirac_trace(const ex & e, unsigned char rl, const ex & trONE)
                        if (num == 5)
                                return trONE * I * eps0123(e.op(1).op(1), e.op(2).op(1), e.op(3).op(1), e.op(4).op(1));
 
-                       // Tr gamma5 gamma.mu1 gamma.mu2 gamma.mu3 gamma.mu4 gamma.mu5 gamma.mu6 = ...
-                       if (num == 7) {
-                               ex i1 = e.op(1).op(1), i2 = e.op(2).op(1),
-                                  i3 = e.op(3).op(1), i4 = e.op(4).op(1),
-                                  i5 = e.op(5).op(1), i6 = e.op(6).op(1);
-                               return trONE * I * (lorentz_g(i1, i2) * eps0123(i3, i4, i5, i6)
-                                                 - lorentz_g(i1, i3) * eps0123(i2, i4, i5, i6)
-                                                 + lorentz_g(i1, i4) * eps0123(i2, i3, i5, i6)
-                                                 - lorentz_g(i1, i5) * eps0123(i2, i3, i4, i6)
-                                                 + lorentz_g(i1, i6) * eps0123(i2, i3, i4, i5)
-                                                 + lorentz_g(i2, i3) * eps0123(i1, i4, i5, i6)
-                                                 - lorentz_g(i2, i4) * eps0123(i1, i3, i5, i6)
-                                                 + lorentz_g(i2, i5) * eps0123(i1, i3, i4, i6)
-                                                 - lorentz_g(i2, i6) * eps0123(i1, i3, i4, i5)
-                                                 + lorentz_g(i3, i4) * eps0123(i1, i2, i5, i6)
-                                                 - lorentz_g(i3, i5) * eps0123(i1, i2, i4, i6)
-                                                 + lorentz_g(i3, i6) * eps0123(i1, i2, i4, i5)
-                                                 + lorentz_g(i4, i5) * eps0123(i1, i2, i3, i6)
-                                                 - lorentz_g(i4, i6) * eps0123(i1, i2, i3, i5)
-                                                 + lorentz_g(i5, i6) * eps0123(i1, i2, i3, i4));
-                       }
-
-                       // Tr gamma5 S_2k =
+                       // Tr gamma5 S_2k =
                        //   I/4! * epsilon0123.mu1.mu2.mu3.mu4 * Tr gamma.mu1 gamma.mu2 gamma.mu3 gamma.mu4 S_2k
+                       exvector ix;
+                       ix.reserve(num - 1);
+                       for (unsigned i=1; i<num; i++)
+                               ix.push_back(e.op(i).op(1));
+                       num--;
+                       int *iv = new int[num];
                        ex result;
-                       for (int i=1; i<num-3; i++) {
-                               ex idx1 = e.op(i).op(1);
+                       for (int i=0; i<num-3; i++) {
+                               ex idx1 = ix[i];
                                for (int j=i+1; j<num-2; j++) {
-                                       ex idx2 = e.op(j).op(1);
+                                       ex idx2 = ix[j];
                                        for (int k=j+1; k<num-1; k++) {
-                                               ex idx3 = e.op(k).op(1);
+                                               ex idx3 = ix[k];
                                                for (int l=k+1; l<num; l++) {
-                                                       ex idx4 = e.op(l).op(1);
-                                                       vector<int> iv;
-                                                       iv.reserve(num-1);
+                                                       ex idx4 = ix[l];
+                                                       iv[0] = i; iv[1] = j; iv[2] = k; iv[3] = l;
                                                        exvector v;
-                                                       v.reserve(num-1);
-                                                       iv.push_back(i); iv.push_back(j); iv.push_back(k); iv.push_back(l);
-                                                       for (int n=1; n<num; n++) {
+                                                       v.reserve(num - 4);
+                                                       for (int n=0, t=4; n<num; n++) {
                                                                if (n == i || n == j || n == k || n == l)
                                                                        continue;
-                                                               iv.push_back(n);
-                                                               v.push_back(e.op(n));
+                                                               iv[t++] = n;
+                                                               v.push_back(ix[n]);
                                                        }
-                                                       int sign = permutation_sign(iv.begin(), iv.end());
+                                                       int sign = permutation_sign(iv, iv + num);
                                                        result += sign * eps0123(idx1, idx2, idx3, idx4)
-                                                               * dirac_trace(ncmul(v, true), rl, trONE);
+                                                               * trace_string(v.begin(), num - 4);
                                                }
                                        }
                                }
                        }
-                       return result * I;
+                       delete[] iv;
+                       return trONE * I * result;
 
                } else { // no gamma5
 
@@ -465,33 +484,12 @@ ex dirac_trace(const ex & e, unsigned char rl, const ex & trONE)
                        if (num == 2)
                                return trONE * lorentz_g(e.op(0).op(1), e.op(1).op(1));
 
-                       // Tr gamma.mu gamma.nu gamma.rho gamma.sig = 4 (g.mu.nu g.rho.sig + g.nu.rho g.mu.sig - g.mu.rho g.nu.sig
-                       if (num == 4)
-                               return trONE * (lorentz_g(e.op(0).op(1), e.op(1).op(1)) * lorentz_g(e.op(2).op(1), e.op(3).op(1))
-                                             + lorentz_g(e.op(1).op(1), e.op(2).op(1)) * lorentz_g(e.op(0).op(1), e.op(3).op(1))
-                                             - lorentz_g(e.op(0).op(1), e.op(2).op(1)) * lorentz_g(e.op(1).op(1), e.op(3).op(1)));
-
-                       // Traces of 6 or more gammas are computed recursively:
-                       // Tr gamma.mu1 gamma.mu2 ... gamma.mun =
-                       //   + g.mu1.mu2 * Tr gamma.mu3 ... gamma.mun
-                       //   - g.mu1.mu3 * Tr gamma.mu2 gamma.mu4 ... gamma.mun
-                       //   + g.mu1.mu4 * Tr gamma.mu3 gamma.mu3 gamma.mu5 ... gamma.mun
-                       //   - ...
-                       //   + g.mu1.mun * Tr gamma.mu2 ... gamma.mu(n-1)
-                       exvector v(num - 2);
-                       int sign = 1;
-                       const ex &ix1 = e.op(0).op(1);
-                       ex result;
-                       for (int i=1; i<num; i++) {
-                               for (int n=1, j=0; n<num; n++) {
-                                       if (n == i)
-                                               continue;
-                                       v[j++] = e.op(n);
-                               }
-                               result += sign * lorentz_g(ix1, e.op(i).op(1)) * dirac_trace(ncmul(v), rl, trONE);
-                               sign = -sign;
-                       }
-                       return result;
+                       exvector iv;
+                       iv.reserve(num);
+                       for (unsigned i=0; i<num; i++)
+                               iv.push_back(e.op(i).op(1));
+
+                       return trONE * trace_string(iv.begin(), num);
                }
        }
 

diff --git a/ginac/indexed.cpp b/ginac/indexed.cpp
index 20e54752a93d89c90e17607c953f0e62741fbacd..f0726915674bf13ce8df0383013ba0bb4ed7fe24 100644 (file)
--- a/ginac/indexed.cpp
+++ b/ginac/indexed.cpp
@@ -308,15 +308,15 @@ ex indexed::eval(int level) const
 
        // If the base object is a product, pull out the numeric factor
        if (is_ex_exactly_of_type(base, mul) && is_ex_exactly_of_type(base.op(base.nops() - 1), numeric)) {
-               exvector v = seq;
+               exvector v(seq);
                ex f = ex_to_numeric(base.op(base.nops() - 1));
                v[0] = seq[0] / f;
                return f * thisexprseq(v);
        }
 
        // Canonicalize indices according to the symmetry properties
-       if (seq.size() > 2 && (symmetry != unknown && symmetry != mixed)) {
-               exvector v = seq;
+       if (seq.size() > 2 && (symmetry == symmetric || symmetry == antisymmetric)) {
+               exvector v(seq);
                int sig = canonicalize_indices(v.begin() + 1, v.end(), symmetry == antisymmetric);
                if (sig != INT_MAX) {
                        // Something has changed while sorting indices, more evaluations later

diff --git a/ginac/inifcns.cpp b/ginac/inifcns.cpp
index 2667215416129be392b21891aa2d16ed6bf89139..beb51d6289dcb56d63e3b12f6d1762ed399e71a9 100644 (file)
--- a/ginac/inifcns.cpp
+++ b/ginac/inifcns.cpp
@@ -533,10 +533,9 @@ static ex symm(const ex & e, exvector::const_iterator first, exvector::const_ite
 
        // 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());
+       iv_lst.insert(iv_lst.begin(), first, last);
+       shaker_sort(iv_lst.begin(), iv_lst.end(), ex_is_less());
        lst orig_lst(iv_lst);
 
        // With n objects there are n! possible permutations

diff --git a/ginac/utils.h b/ginac/utils.h
index efdeaa2b7fe5cee59435176dcb32d216286f27c9..c0bc2efecf49c00ec10328bc5b51674b00341acb 100644 (file)
--- a/ginac/utils.h
+++ b/ginac/utils.h
@@ -128,68 +128,111 @@ inline unsigned golden_ratio_hash(unsigned n)
 }
 
 /* Compute the sign of a permutation of a container, with and without an
-   explicitly supplied comparison function. The containers gets modified
-   during the operation. */
+   explicitly supplied comparison function. If the sign returned is 1 or -1,
+   the container is sorted after the operation. */
 template <class It>
-int permutation_sign(It first, It last)
+inline int permutation_sign(It first, It last)
 {
        if (first == last)
                return 0;
-       It i = first;
-       ++i;
-       if (i == last)
+       --last;
+       if (first == last)
                return 0;
-       i = first;
-       It next_to_last = last;
-       --next_to_last;
-
+       It flag = first;
        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);
+
+       do {
+               It i = last, other = last;
+               --other;
+               bool swapped = false;
+               while (i != first) {
+                       if (*i < *other) {
+                               iter_swap(other, i);
+                               flag = other;
+                               swapped = true;
                                sign = -sign;
-                       }
-                       ++j;
+                       } else if (!(*other < *i))
+                               return 0;
+                       --i; --other;
+               }
+               if (!swapped)
+                       return sign;
+               ++flag;
+               if (flag == last)
+                       return sign;
+               first = flag;
+               i = first; other = first;
+               ++other;
+               swapped = false;
+               while (i != last) {
+                       if (*other < *i) {
+                               iter_swap(i, other);
+                               flag = other;
+                               swapped = true;
+                               sign = -sign;
+                       } else if (!(*i < *other))
+                               return 0;
+                       ++i; ++other;
                }
-               ++i;
-       }
+               if (!swapped)
+                       return sign;
+               last = flag;
+               --last;
+       } while (first != last);
+
        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)
+inline int permutation_sign(It first, It last, Cmp comp)
 {
        if (first == last)
                return 0;
-       It i = first;
-       ++i;
-       if (i == last)
+       --last;
+       if (first == last)
                return 0;
-       i = first;
-       It next_to_last = last;
-       --next_to_last;
-
+       It flag = first;
        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);
+
+       do {
+               It i = last, other = last;
+               --other;
+               bool swapped = false;
+               while (i != first) {
+                       if (comp(*i, *other)) {
+                               iter_swap(other, i);
+                               flag = other;
+                               swapped = true;
                                sign = -sign;
-                       }
-                       ++j;
+                       } else if (!comp(*other, *i))
+                               return 0;
+                       --i; --other;
                }
-               ++i;
-       }
+               if (!swapped)
+                       return sign;
+               ++flag;
+               if (flag == last)
+                       return sign;
+               first = flag;
+               i = first; other = first;
+               ++other;
+               swapped = false;
+               while (i != last) {
+                       if (comp(*other, *i)) {
+                               iter_swap(i, other);
+                               flag = other;
+                               swapped = true;
+                               sign = -sign;
+                       } else if (!comp(*i, *other))
+                               return 0;
+                       ++i; ++other;
+               }
+               if (!swapped)
+                       return sign;
+               last = flag;
+               --last;
+       } while (first != last);
+
        return sign;
 }
 
@@ -203,30 +246,41 @@ inline void shaker_sort(It first, It last, Cmp comp)
        if (first == last)
                return;
        It flag = first;
+
        do {
                It i = last, other = last;
                --other;
-               while (i > first) {
+               bool swapped = false;
+               while (i != first) {
                        if (comp(*i, *other)) {
                                iter_swap(other, i);
                                flag = other;
+                               swapped = true;
                        }
                        --i; --other;
                }
+               if (!swapped)
+                       return;
                ++flag;
+               if (flag == last)
+                       return;
                first = flag;
                i = first; other = first;
                ++other;
-               while (i < last) {
+               swapped = false;
+               while (i != last) {
                        if (comp(*other, *i)) {
                                iter_swap(i, other);
                                flag = other;
+                               swapped = true;
                        }
                        ++i; ++other;
                }
+               if (!swapped)
+                       return;
                last = flag;
                --last;
-       } while (first <= last);
+       } while (first != last);
 }
 
 /* Function objects for STL sort() etc. */