- added variants of dirac_trace() and color_trace() that take the trace over

  more than one representation label by specifying a set or list of labels
- diracgamma::contract_with() uses Chisholm identities in 4 dimensions to
  produce more compact results
- added missing documentation for options parameter of simplify_indexed()

diff --git a/ginac/clifford.cpp b/ginac/clifford.cpp
index f8d503ff6e68ea15f8bdcd8fdc107d85e72f2081..b2e56224e11ee9518ccfe87fdffea1968686591c 100644 (file)
--- a/ginac/clifford.cpp
+++ b/ginac/clifford.cpp
@@ -250,6 +250,14 @@ static void base_and_index(const ex & c, ex & b, ex & i)
        }
 }
 
        }
 }
 

+/** Predicate for finding non-clifford objects. */
+struct is_not_a_clifford : public std::unary_function<ex, bool> {
+       bool operator()(const ex & e)
+       {
+               return !is_a<clifford>(e);
+       }
+};
+

 /** Contraction of a gamma matrix with something else. */
 bool diracgamma::contract_with(exvector::iterator self, exvector::iterator other, exvector & v) const
 {
 /** Contraction of a gamma matrix with something else. */
 bool diracgamma::contract_with(exvector::iterator self, exvector::iterator other, exvector & v) const
 {


@@ -268,21 +276,23 @@ bool diracgamma::contract_with(exvector::iterator self, exvector::iterator other
                if (ex_to<clifford>(*other).get_representation_label() != rl)
                        return false;
 
                if (ex_to<clifford>(*other).get_representation_label() != rl)
                        return false;
 

+               size_t num = other - self;
+

                // gamma~mu gamma.mu = dim ONE
                // gamma~mu gamma.mu = dim ONE

-               if (other - self == 1) {
+               if (num == 1) {

                        *self = dim;
                        *other = dirac_ONE(rl);
                        return true;
 
                // gamma~mu gamma~alpha gamma.mu = (2-dim) gamma~alpha
                        *self = dim;
                        *other = dirac_ONE(rl);
                        return true;
 
                // gamma~mu gamma~alpha gamma.mu = (2-dim) gamma~alpha

-               } else if (other - self == 2
+               } else if (num == 2

                        && is_a<clifford>(self[1])) {
                        *self = 2 - dim;
                        *other = _ex1;
                        return true;
 
                // gamma~mu gamma~alpha gamma~beta gamma.mu = 4 g~alpha~beta + (dim-4) gamam~alpha gamma~beta
                        && is_a<clifford>(self[1])) {
                        *self = 2 - dim;
                        *other = _ex1;
                        return true;
 
                // gamma~mu gamma~alpha gamma~beta gamma.mu = 4 g~alpha~beta + (dim-4) gamam~alpha gamma~beta

-               } else if (other - self == 3
+               } else if (num == 3

                        && is_a<clifford>(self[1])
                        && is_a<clifford>(self[2])) {
                        ex b1, i1, b2, i2;
                        && is_a<clifford>(self[1])
                        && is_a<clifford>(self[2])) {
                        ex b1, i1, b2, i2;


@@ -295,7 +305,7 @@ bool diracgamma::contract_with(exvector::iterator self, exvector::iterator other
                        return true;
 
                // gamma~mu gamma~alpha gamma~beta gamma~delta gamma.mu = -2 gamma~delta gamma~beta gamma~alpha - (dim-4) gamam~alpha gamma~beta gamma~delta
                        return true;
 
                // gamma~mu gamma~alpha gamma~beta gamma~delta gamma.mu = -2 gamma~delta gamma~beta gamma~alpha - (dim-4) gamam~alpha gamma~beta gamma~delta

-               } else if (other - self == 4
+               } else if (num == 4

                        && is_a<clifford>(self[1])
                        && is_a<clifford>(self[2])
                        && is_a<clifford>(self[3])) {
                        && is_a<clifford>(self[1])
                        && is_a<clifford>(self[2])
                        && is_a<clifford>(self[3])) {


@@ -306,27 +316,45 @@ bool diracgamma::contract_with(exvector::iterator self, exvector::iterator other
                        *other = _ex1;
                        return true;
 
                        *other = _ex1;
                        return true;
 

+               // gamma~mu Sodd gamma.mu = -2 Sodd_R
+               // (Chisholm identity in 4 dimensions)
+               } else if (!((other - self) & 1) && dim.is_equal(4)) {
+                       if (std::find_if(self + 1, other, is_not_a_clifford()) != other)
+                               return false;
+
+                       *self = ncmul(exvector(std::reverse_iterator<exvector::const_iterator>(other), std::reverse_iterator<exvector::const_iterator>(self + 1)), true);
+                       std::fill(self + 1, other, _ex1);
+                       *other = _ex_2;
+                       return true;
+
+               // gamma~mu Sodd gamma~alpha gamma.mu = 2 gamma~alpha Sodd + 2 Sodd_R gamma~alpha
+               // (commutate contracted indices towards each other, then use
+               // Chisholm identity in 4 dimensions)
+               } else if (((other - self) & 1) && dim.is_equal(4)) {
+                       if (std::find_if(self + 1, other, is_not_a_clifford()) != other)
+                               return false;
+
+                       exvector::iterator next_to_last = other - 1;
+                       ex S = ncmul(exvector(self + 1, next_to_last), true);
+                       ex SR = ncmul(exvector(std::reverse_iterator<exvector::const_iterator>(next_to_last), std::reverse_iterator<exvector::const_iterator>(self + 1)), true);
+
+                       *self = (*next_to_last) * S + SR * (*next_to_last);
+                       std::fill(self + 1, other, _ex1);
+                       *other = _ex2;
+                       return true;
+

                // gamma~mu S gamma~alpha gamma.mu = 2 gamma~alpha S - gamma~mu S gamma.mu gamma~alpha
                // (commutate contracted indices towards each other, simplify_indexed()
                // will re-expand and re-run the simplification)
                } else {
                // gamma~mu S gamma~alpha gamma.mu = 2 gamma~alpha S - gamma~mu S gamma.mu gamma~alpha
                // (commutate contracted indices towards each other, simplify_indexed()
                // will re-expand and re-run the simplification)
                } else {

-                       exvector::iterator it = self + 1, next_to_last = other - 1;
-                       while (it != other) {
-                               if (!is_a<clifford>(*it))
-                                       return false;
-                               ++it;
-                       }
+                       if (std::find_if(self + 1, other, is_not_a_clifford()) != other)
+                               return false;

 
 

-                       it = self + 1;
-                       ex S = _ex1;
-                       while (it != next_to_last) {
-                               S *= *it;
-                               *it++ = _ex1;
-                       }
+                       exvector::iterator next_to_last = other - 1;
+                       ex S = ncmul(exvector(self + 1, next_to_last), true);

 
                        *self = 2 * (*next_to_last) * S - (*self) * S * (*other) * (*next_to_last);
 
                        *self = 2 * (*next_to_last) * S - (*self) * S * (*other) * (*next_to_last);

-                       *next_to_last = _ex1;
-                       *other = _ex1;
+                       std::fill(self + 1, other + 1, _ex1);

                        return true;
                }
 
                        return true;
                }
 


@@ -702,6 +730,13 @@ static bool is_clifford_tinfo(unsigned ti)
        return (ti & ~0xff) == TINFO_clifford;
 }
 
        return (ti & ~0xff) == TINFO_clifford;
 }
 

+/** Extract representation label from tinfo key (as returned by
+ *  return_type_tinfo()). */
+static unsigned char get_representation_label(unsigned ti)
+{
+       return ti & 0xff;
+}
+

 /** 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, size_t num)
 /** 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, size_t num)


@@ -738,12 +773,17 @@ static ex trace_string(exvector::const_iterator ix, size_t num)
        return result;
 }
 
        return result;
 }
 

-ex dirac_trace(const ex & e, unsigned char rl, const ex & trONE)
+ex dirac_trace(const ex & e, const std::set<unsigned char> & rls, const ex & trONE)

 {
        if (is_a<clifford>(e)) {
 
 {
        if (is_a<clifford>(e)) {
 

-               if (!ex_to<clifford>(e).get_representation_label() == rl)
-                       return _ex0;
+               unsigned char rl = ex_to<clifford>(e).get_representation_label();
+
+               // Are we taking the trace over this object's representation label?
+               if (rls.find(rl) == rls.end())
+                       return e;
+
+               // Yes, all elements are traceless, except for dirac_ONE and dirac_L/R

                const ex & g = e.op(0);
                if (is_a<diracone>(g))
                        return trONE;
                const ex & g = e.op(0);
                if (is_a<diracone>(g))
                        return trONE;


@@ -758,8 +798,8 @@ ex dirac_trace(const ex & e, unsigned char rl, const ex & trONE)
                ex prod = _ex1;
                for (size_t i=0; i<e.nops(); i++) {
                        const ex &o = e.op(i);
                ex prod = _ex1;
                for (size_t i=0; i<e.nops(); i++) {
                        const ex &o = e.op(i);

-                       if (is_clifford_tinfo(o.return_type_tinfo(), rl))
-                               prod *= dirac_trace(o, rl, trONE);
+                       if (is_clifford_tinfo(o.return_type_tinfo()))
+                               prod *= dirac_trace(o, rls, trONE);

                        else
                                prod *= o;
                }
                        else
                                prod *= o;
                }


@@ -767,8 +807,11 @@ ex dirac_trace(const ex & e, unsigned char rl, const ex & trONE)
 
        } else if (is_exactly_a<ncmul>(e)) {
 
 
        } else if (is_exactly_a<ncmul>(e)) {
 

-               if (!is_clifford_tinfo(e.return_type_tinfo(), rl))
-                       return _ex0;
+               unsigned char rl = get_representation_label(e.return_type_tinfo());
+
+               // Are we taking the trace over this string's representation label?
+               if (rls.find(rl) == rls.end())
+                       return e;

 
                // Substitute gammaL/R and expand product, if necessary
                ex e_expanded = e.subs(lst(
 
                // Substitute gammaL/R and expand product, if necessary
                ex e_expanded = e.subs(lst(


@@ -776,7 +819,7 @@ ex dirac_trace(const ex & e, unsigned char rl, const ex & trONE)
                        dirac_gammaR(rl) == (dirac_ONE(rl)+dirac_gamma5(rl))/2
                ), subs_options::no_pattern).expand();
                if (!is_a<ncmul>(e_expanded))
                        dirac_gammaR(rl) == (dirac_ONE(rl)+dirac_gamma5(rl))/2
                ), subs_options::no_pattern).expand();
                if (!is_a<ncmul>(e_expanded))

-                       return dirac_trace(e_expanded, rl, trONE);
+                       return dirac_trace(e_expanded, rls, trONE);

 
                // gamma5 gets moved to the front so this check is enough
                bool has_gamma5 = is_a<diracgamma5>(e.op(0).op(0));
 
                // gamma5 gets moved to the front so this check is enough
                bool has_gamma5 = is_a<diracgamma5>(e.op(0).op(0));


@@ -860,13 +903,35 @@ ex dirac_trace(const ex & e, unsigned char rl, const ex & trONE)
        } else if (e.nops() > 0) {
 
                // Trace maps to all other container classes (this includes sums)
        } else if (e.nops() > 0) {
 
                // Trace maps to all other container classes (this includes sums)

-               pointer_to_map_function_2args<unsigned char, const ex &> fcn(dirac_trace, rl, trONE);
+               pointer_to_map_function_2args<const std::set<unsigned char> &, const ex &> fcn(dirac_trace, rls, trONE);

                return e.map(fcn);
 
        } else
                return _ex0;
 }
 
                return e.map(fcn);
 
        } else
                return _ex0;
 }
 

+ex dirac_trace(const ex & e, const lst & rll, const ex & trONE)
+{
+       // Convert list to set
+       std::set<unsigned char> rls;
+       for (lst::const_iterator i = rll.begin(); i != rll.end(); ++i) {
+               if (i->info(info_flags::nonnegint))
+                       rls.insert(ex_to<numeric>(*i).to_int());
+       }
+
+       return dirac_trace(e, rls, trONE);
+}
+
+ex dirac_trace(const ex & e, unsigned char rl, const ex & trONE)
+{
+       // Convert label to set
+       std::set<unsigned char> rls;
+       rls.insert(rl);
+
+       return dirac_trace(e, rls, trONE);
+}
+
+

 ex canonicalize_clifford(const ex & e)
 {
        // Scan for any ncmul objects
 ex canonicalize_clifford(const ex & e)
 {
        // Scan for any ncmul objects

diff --git a/ginac/clifford.h b/ginac/clifford.h
index 667d7befd9d400749188f8af6f28bcecee1ce83e..81abd9f5c879c956da862f009fc5f24f6b7bcbc4 100644 (file)
--- a/ginac/clifford.h
+++ b/ginac/clifford.h
@@ -28,6 +28,8 @@
 #include "symbol.h"
 #include "idx.h"
 
 #include "symbol.h"
 #include "idx.h"
 

+#include <set>
+

 namespace GiNaC {
 
 
 namespace GiNaC {
 
 


@@ -226,6 +228,26 @@ ex dirac_gammaR(unsigned char rl = 0);
  *  @param rl Representation label */
 ex dirac_slash(const ex & e, const ex & dim, unsigned char rl = 0);
 
  *  @param rl Representation label */
 ex dirac_slash(const ex & e, const ex & dim, unsigned char rl = 0);
 

+/** Calculate dirac traces over the specified set of representation labels.
+ *  The computed trace is a linear functional that is equal to the usual
+ *  trace only in D = 4 dimensions. In particular, the functional is not
+ *  always cyclic in D != 4 dimensions when gamma5 is involved.
+ *
+ *  @param e Expression to take the trace of
+ *  @param rls Set of representation labels
+ *  @param trONE Expression to be returned as the trace of the unit matrix */
+ex dirac_trace(const ex & e, const std::set<unsigned char> & rls, const ex & trONE = 4);
+
+/** Calculate dirac traces over the specified list of representation labels.
+ *  The computed trace is a linear functional that is equal to the usual
+ *  trace only in D = 4 dimensions. In particular, the functional is not
+ *  always cyclic in D != 4 dimensions when gamma5 is involved.
+ *
+ *  @param e Expression to take the trace of
+ *  @param rll List of representation labels
+ *  @param trONE Expression to be returned as the trace of the unit matrix */
+ex dirac_trace(const ex & e, const lst & rll, const ex & trONE = 4);
+

 /** Calculate the trace of an expression containing gamma objects with
  *  a specified representation label. The computed trace is a linear
  *  functional that is equal to the usual trace only in D = 4 dimensions.
 /** Calculate the trace of an expression containing gamma objects with
  *  a specified representation label. The computed trace is a linear
  *  functional that is equal to the usual trace only in D = 4 dimensions.

diff --git a/ginac/color.cpp b/ginac/color.cpp
index e06eff44a43c659ef7bb7567f6dc04fdffc03961..3f45bba5acaaa4bf9983245f650442f46091e3a2 100644 (file)
--- a/ginac/color.cpp
+++ b/ginac/color.cpp
@@ -527,12 +527,32 @@ static bool is_color_tinfo(unsigned ti, unsigned char rl)
        return ti == (TINFO_color + rl);
 }
 
        return ti == (TINFO_color + rl);
 }
 

-ex color_trace(const ex & e, unsigned char rl)
+/** Check whether a given tinfo key (as returned by return_type_tinfo()
+ *  is that of a color object (with an arbitrary representation label). */
+static bool is_color_tinfo(unsigned ti)
+{
+       return (ti & ~0xff) == TINFO_color;
+}
+
+/** Extract representation label from tinfo key (as returned by
+ *  return_type_tinfo()). */
+static unsigned char get_representation_label(unsigned ti)
+{
+       return ti & 0xff;
+}
+
+ex color_trace(const ex & e, const std::set<unsigned char> & rls)

 {
        if (is_a<color>(e)) {
 
 {
        if (is_a<color>(e)) {
 

-               if (ex_to<color>(e).get_representation_label() == rl
-                && is_a<su3one>(e.op(0)))
+               unsigned char rl = ex_to<color>(e).get_representation_label();
+
+               // Are we taking the trace over this object's representation label?
+               if (rls.find(rl) == rls.end())
+                       return e;
+
+               // Yes, all generators are traceless, except for color_ONE
+               if (is_a<su3one>(e.op(0)))

                        return _ex3;
                else
                        return _ex0;
                        return _ex3;
                else
                        return _ex0;


@@ -543,8 +563,8 @@ ex color_trace(const ex & e, unsigned char rl)
                ex prod = _ex1;
                for (size_t i=0; i<e.nops(); i++) {
                        const ex &o = e.op(i);
                ex prod = _ex1;
                for (size_t i=0; i<e.nops(); i++) {
                        const ex &o = e.op(i);

-                       if (is_color_tinfo(o.return_type_tinfo(), rl))
-                               prod *= color_trace(o, rl);
+                       if (is_color_tinfo(o.return_type_tinfo()))
+                               prod *= color_trace(o, rls);

                        else
                                prod *= o;
                }
                        else
                                prod *= o;
                }


@@ -552,13 +572,16 @@ ex color_trace(const ex & e, unsigned char rl)
 
        } else if (is_exactly_a<ncmul>(e)) {
 
 
        } else if (is_exactly_a<ncmul>(e)) {
 

-               if (!is_color_tinfo(e.return_type_tinfo(), rl))
-                       return _ex0;
+               unsigned char rl = get_representation_label(e.return_type_tinfo());

 
 

-               // Expand product, if necessary
+               // Are we taking the trace over this string's representation label?
+               if (rls.find(rl) == rls.end())
+                       return e;
+
+               // Yes, expand product if necessary

                ex e_expanded = e.expand();
                if (!is_a<ncmul>(e_expanded))
                ex e_expanded = e.expand();
                if (!is_a<ncmul>(e_expanded))

-                       return color_trace(e_expanded, rl);
+                       return color_trace(e_expanded, rls);

 
                size_t num = e.nops();
 
 
                size_t num = e.nops();
 


@@ -597,11 +620,32 @@ ex color_trace(const ex & e, unsigned char rl)
        } else if (e.nops() > 0) {
 
                // Trace maps to all other container classes (this includes sums)
        } else if (e.nops() > 0) {
 
                // Trace maps to all other container classes (this includes sums)

-               pointer_to_map_function_1arg<unsigned char> fcn(color_trace, rl);
+               pointer_to_map_function_1arg<const std::set<unsigned char> &> fcn(color_trace, rls);

                return e.map(fcn);
 
        } else
                return _ex0;
 }
 
                return e.map(fcn);
 
        } else
                return _ex0;
 }
 

+ex color_trace(const ex & e, const lst & rll)
+{
+       // Convert list to set
+       std::set<unsigned char> rls;
+       for (lst::const_iterator i = rll.begin(); i != rll.end(); ++i) {
+               if (i->info(info_flags::nonnegint))
+                       rls.insert(ex_to<numeric>(*i).to_int());
+       }
+
+       return color_trace(e, rls);
+}
+
+ex color_trace(const ex & e, unsigned char rl)
+{
+       // Convert label to set
+       std::set<unsigned char> rls;
+       rls.insert(rl);
+
+       return color_trace(e, rls);
+}
+

 } // namespace GiNaC
 } // namespace GiNaC

diff --git a/ginac/color.h b/ginac/color.h
index 2086909e1304dd1c7a35ac9e53c658263150348c..826b4c7ffe68ae7489febb13aa82db7e705ab0f2 100644 (file)
--- a/ginac/color.h
+++ b/ginac/color.h
@@ -26,6 +26,8 @@
 #include "indexed.h"
 #include "tensor.h"
 
 #include "indexed.h"
 #include "tensor.h"
 

+#include <set>
+

 namespace GiNaC {
 
 
 namespace GiNaC {
 
 


@@ -169,6 +171,18 @@ ex color_d(const ex & a, const ex & b, const ex & c);
 /** This returns the linear combination d.a.b.c+I*f.a.b.c. */
 ex color_h(const ex & a, const ex & b, const ex & c);
 
 /** This returns the linear combination d.a.b.c+I*f.a.b.c. */
 ex color_h(const ex & a, const ex & b, const ex & c);
 

+/** Calculate color traces over the specified set of representation labels.
+ *
+ *  @param e Expression to take the trace of
+ *  @param rls Set of representation labels */
+ex color_trace(const ex & e, const std::set<unsigned char> & rls);
+
+/** Calculate color traces over the specified list of representation labels.
+ *
+ *  @param e Expression to take the trace of
+ *  @param rll List of representation labels */
+ex color_trace(const ex & e, const lst & rll);
+

 /** Calculate the trace of an expression containing color objects with a
  *  specified representation label.
  *
 /** Calculate the trace of an expression containing color objects with a
  *  specified representation label.
  *

diff --git a/ginac/indexed.cpp b/ginac/indexed.cpp
index ad17a497a9dc703e52cd435d6901312071beab42..d693373cae62c6cd622495d05722aa9ce7065d99 100644 (file)
--- a/ginac/indexed.cpp
+++ b/ginac/indexed.cpp
@@ -1128,6 +1128,7 @@ ex simplify_indexed(const ex & e, exvector & free_indices, exvector & dummy_indi
  *  performs contraction of dummy indices where possible and checks whether
  *  the free indices in sums are consistent.
  *
  *  performs contraction of dummy indices where possible and checks whether
  *  the free indices in sums are consistent.
  *

+ *  @param options Simplification options (currently unused)

  *  @return simplified expression */
 ex ex::simplify_indexed(unsigned options) const
 {
  *  @return simplified expression */
 ex ex::simplify_indexed(unsigned options) const
 {


@@ -1142,6 +1143,7 @@ ex ex::simplify_indexed(unsigned options) const
  *  scalar products by known values if desired.
  *
  *  @param sp Scalar products to be replaced automatically
  *  scalar products by known values if desired.
  *
  *  @param sp Scalar products to be replaced automatically

+ *  @param options Simplification options (currently unused)

  *  @return simplified expression */
 ex ex::simplify_indexed(const scalar_products & sp, unsigned options) const
 {
  *  @return simplified expression */
 ex ex::simplify_indexed(const scalar_products & sp, unsigned options) const
 {