- added Clifford algebra unity element

- superfluous unity elements are removed from Clifford and color strings
- added Clifford contractions:
    gamma~mu*gamma.mu
    gamma~mu*gamma~alpha*gamma.mu
- added color contractions:
    d.aac
    d.abc*d.abc
    d.akl*d.bkl
    d.abc*f.abc
    d.akl*f.bkl
    f.abc*f.abc
    f.akl*f.bkl
- delta tensor and color structure constants can be evaluated numerically
- color_T(), color_d() and color_f() check their arguments
- added a couple of utility functions to idx.*
- simplify_indexed() doesn't crash any more when used on expressions containing
  noncommutative products

diff --git a/ginac/clifford.cpp b/ginac/clifford.cpp
index 5875bbb836e80a7c4ed80789b5c934c0d618f469..94f434a7710da07db1625332d7fae96b8968c095 100644 (file)
--- a/ginac/clifford.cpp
+++ b/ginac/clifford.cpp
@@ -33,6 +33,7 @@
 namespace GiNaC {
 
 GINAC_IMPLEMENT_REGISTERED_CLASS(clifford, indexed)
 namespace GiNaC {
 
 GINAC_IMPLEMENT_REGISTERED_CLASS(clifford, indexed)

+GINAC_IMPLEMENT_REGISTERED_CLASS(diracone, tensor)

 GINAC_IMPLEMENT_REGISTERED_CLASS(diracgamma, tensor)
 
 //////////
 GINAC_IMPLEMENT_REGISTERED_CLASS(diracgamma, tensor)
 
 //////////


@@ -47,6 +48,7 @@ clifford::clifford()
 
 DEFAULT_COPY(clifford)
 DEFAULT_DESTROY(clifford)
 
 DEFAULT_COPY(clifford)
 DEFAULT_DESTROY(clifford)

+DEFAULT_CTORS(diracone)

 DEFAULT_CTORS(diracgamma)
 
 //////////
 DEFAULT_CTORS(diracgamma)
 
 //////////


@@ -63,6 +65,15 @@ clifford::clifford(const ex & b, const ex & mu) : inherited(b, mu)
        tinfo_key = TINFO_clifford;
 }
 
        tinfo_key = TINFO_clifford;
 }
 

+/** Construct object without any indices. This constructor is for internal
+ *  use only. Use the dirac_one() function instead.
+ *  @see dirac_one */
+clifford::clifford(const ex & b) : inherited(b)
+{
+       debugmsg("clifford constructor from ex", LOGLEVEL_CONSTRUCT);
+       tinfo_key = TINFO_clifford;
+}
+

 clifford::clifford(const exvector & v, bool discardable) : inherited(indexed::unknown, v, discardable)
 {
        debugmsg("clifford constructor from exvector", LOGLEVEL_CONSTRUCT);
 clifford::clifford(const exvector & v, bool discardable) : inherited(indexed::unknown, v, discardable)
 {
        debugmsg("clifford constructor from exvector", LOGLEVEL_CONSTRUCT);


@@ -80,6 +91,7 @@ clifford::clifford(exvector * vp) : inherited(indexed::unknown, vp)
 //////////
 
 DEFAULT_ARCHIVING(clifford)
 //////////
 
 DEFAULT_ARCHIVING(clifford)

+DEFAULT_ARCHIVING(diracone)

 DEFAULT_ARCHIVING(diracgamma)
 
 //////////
 DEFAULT_ARCHIVING(diracgamma)
 
 //////////


@@ -91,15 +103,57 @@ int clifford::compare_same_type(const basic & other) const
        return inherited::compare_same_type(other);
 }
 
        return inherited::compare_same_type(other);
 }
 

+DEFAULT_COMPARE(diracone)

 DEFAULT_COMPARE(diracgamma)
 DEFAULT_COMPARE(diracgamma)

+DEFAULT_PRINT(diracone, "ONE")

 DEFAULT_PRINT(diracgamma, "gamma")
 
 DEFAULT_PRINT(diracgamma, "gamma")
 

+/** Contraction of a gamma matrix with something else. */
+bool diracgamma::contract_with(exvector::iterator self, exvector::iterator other, exvector & v) const
+{
+       GINAC_ASSERT(is_ex_of_type(*self, indexed));
+       GINAC_ASSERT(is_ex_of_type(*other, indexed));
+       GINAC_ASSERT(is_ex_of_type(self->op(0), diracgamma));
+
+       if (is_ex_of_type(other->op(0), diracgamma)) {
+
+               // gamma~mu*gamma.mu = dim*ONE
+               if (other - self == 1) {
+                       *self = ex_to_idx(self->op(1)).get_dim();
+                       *other = dirac_one();
+                       return true;
+
+               // gamma~mu*gamma~alpha*gamma.mu = (2-dim)*gamma~alpha
+               } else if (other - self == 2) {
+                       *self = 2 - ex_to_idx(self->op(1)).get_dim();
+                       *other = _ex1();
+                       return true;
+               }
+       }
+
+       return false;
+}
+

 /** Perform automatic simplification on noncommutative product of clifford
 /** Perform automatic simplification on noncommutative product of clifford

- *  objects. */
+ *  objects. This removes superfluous ONEs. */

 ex clifford::simplify_ncmul(const exvector & v) const
 {
 ex clifford::simplify_ncmul(const exvector & v) const
 {

-       //!! to be implemented
-       return nonsimplified_ncmul(v);
+       exvector s;
+       s.reserve(v.size());
+
+       exvector::const_iterator it = v.begin(), itend = v.end();
+       while (it != itend) {
+               if (!is_ex_of_type(it->op(0), diracone))
+                       s.push_back(*it);
+               it++;
+       }
+
+       if (s.size() == 0)
+               return clifford(diracone());
+       else if (s.size() == v.size())
+               return simplified_ncmul(v);
+       else
+               return simplified_ncmul(s);

 }
 
 ex clifford::thisexprseq(const exvector & v) const
 }
 
 ex clifford::thisexprseq(const exvector & v) const


@@ -116,10 +170,16 @@ ex clifford::thisexprseq(exvector * vp) const
 // global functions
 //////////
 
 // global functions
 //////////
 

+ex dirac_one(void)
+{
+       return clifford(diracone());
+}
+

 ex dirac_gamma(const ex & mu)
 {
        if (!is_ex_of_type(mu, varidx))
                throw(std::invalid_argument("index of Dirac gamma must be of type varidx"));
 ex dirac_gamma(const ex & mu)
 {
        if (!is_ex_of_type(mu, varidx))
                throw(std::invalid_argument("index of Dirac gamma must be of type varidx"));

+

        return clifford(diracgamma(), mu);
 }
 
        return clifford(diracgamma(), mu);
 }
 

diff --git a/ginac/clifford.h b/ginac/clifford.h
index ccb16bf24c78fa2733a4d759c831a6b1de6af2bc..0feef3f049b135c3984af3bfba8f52df5df70923 100644 (file)
--- a/ginac/clifford.h
+++ b/ginac/clifford.h
@@ -38,6 +38,7 @@ class clifford : public indexed
 
        // other constructors
 public:
 
        // other constructors
 public:

+       clifford(const ex & b);

        clifford(const ex & b, const ex & mu);
 
        // internal constructors
        clifford(const ex & b, const ex & mu);
 
        // internal constructors


@@ -54,6 +55,17 @@ protected:
 };
 
 
 };
 
 

+/** This class represents the Clifford algebra unity element. */
+class diracone : public tensor
+{
+       GINAC_DECLARE_REGISTERED_CLASS(diracone, tensor)
+
+       // functions overriding virtual functions from bases classes
+public:
+       void print(std::ostream & os, unsigned upper_precedence=0) const;
+};
+
+

 /** This class represents the Dirac gamma Lorentz vector. */
 class diracgamma : public tensor
 {
 /** This class represents the Dirac gamma Lorentz vector. */
 class diracgamma : public tensor
 {


@@ -62,6 +74,7 @@ class diracgamma : public tensor
        // functions overriding virtual functions from bases classes
 public:
        void print(std::ostream & os, unsigned upper_precedence=0) const;
        // functions overriding virtual functions from bases classes
 public:
        void print(std::ostream & os, unsigned upper_precedence=0) const;

+       bool contract_with(exvector::iterator self, exvector::iterator other, exvector & v) const;

 };
 
 
 };
 
 


@@ -72,6 +85,11 @@ inline const clifford &ex_to_clifford(const ex &e)
 }
 
 
 }
 
 

+/** Create a Clifford unity object.
+ *
+ *  @return newly constructed object */
+ex dirac_one(void);
+

 /** Create a Dirac gamma object.
  *
  *  @param mu Index (must be of class varidx or a derived class)
 /** Create a Dirac gamma object.
  *
  *  @param mu Index (must be of class varidx or a derived class)

diff --git a/ginac/color.cpp b/ginac/color.cpp
index cd3179ae65cab408a506668859e5b25176bd0898..65dcf4ae436d2648ebf7f83a143c606c90e7e22a 100644 (file)
--- a/ginac/color.cpp
+++ b/ginac/color.cpp
@@ -20,10 +20,15 @@
  *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
  */
 
  *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
  */
 

+#include <algorithm>
+#include <stdexcept>
+

 #include "color.h"
 #include "ex.h"
 #include "color.h"
 #include "ex.h"

+#include "idx.h"

 #include "ncmul.h"
 #include "numeric.h"
 #include "ncmul.h"
 #include "numeric.h"

+#include "power.h" // for sqrt()

 #include "archive.h"
 #include "debugmsg.h"
 #include "utils.h"
 #include "archive.h"
 #include "debugmsg.h"
 #include "utils.h"


@@ -145,8 +150,23 @@ DEFAULT_PRINT(su3d, "d")
  *  objects. This removes superfluous ONEs. */
 ex color::simplify_ncmul(const exvector & v) const
 {
  *  objects. This removes superfluous ONEs. */
 ex color::simplify_ncmul(const exvector & v) const
 {

-       //!! to be implemented
-       return nonsimplified_ncmul(v);
+       //!! TODO: sort by representation label
+       exvector s;
+       s.reserve(v.size());
+
+       exvector::const_iterator it = v.begin(), itend = v.end();
+       while (it != itend) {
+               if (!is_ex_of_type(it->op(0), su3one))
+                       s.push_back(*it);
+               it++;
+       }
+
+       if (s.size() == 0)
+               return color(su3one());
+       else if (s.size() == v.size())
+               return simplified_ncmul(v);
+       else
+               return simplified_ncmul(s);

 }
 
 ex color::thisexprseq(const exvector & v) const
 }
 
 ex color::thisexprseq(const exvector & v) const


@@ -159,6 +179,201 @@ ex color::thisexprseq(exvector * vp) const
        return color(representation_label, vp);
 }
 
        return color(representation_label, vp);
 }
 

+/** Given a vector iv3 of three indices and a vector iv2 of two indices that
+ *  is a subset of iv3, return the (free) index that is in iv3 but not in
+ *  iv2 and the sign introduced by permuting that index to the front.
+ *
+ *  @param iv3 Vector of 3 indices
+ *  @param iv2 Vector of 2 indices, must be a subset of iv3
+ *  @param sig Returs sign introduced by index permutation
+ *  @return the free index (the one that is in iv3 but not in iv2) */
+static ex permute_free_index_to_front(const exvector & iv3, const exvector & iv2, int & sig)
+{
+       GINAC_ASSERT(iv3.size() == 3);
+       GINAC_ASSERT(iv2.size() == 2);
+
+       sig = 1;
+
+#define TEST_PERMUTATION(A,B,C,P) \
+       if (iv3[B].is_equal(iv2[0]) && iv3[C].is_equal(iv2[1])) { \
+               sig = P; \
+               return iv3[A]; \
+       }
+       
+       TEST_PERMUTATION(0,1,2,  1);
+       TEST_PERMUTATION(0,2,1, -1);
+       TEST_PERMUTATION(1,0,2, -1);
+       TEST_PERMUTATION(1,2,0,  1);
+       TEST_PERMUTATION(2,0,1,  1);
+       TEST_PERMUTATION(2,1,0, -1);
+
+       throw(std::logic_error("permute_free_index_to_front(): no valid permutation found"));
+}
+
+/** Automatic symbolic evaluation of indexed symmetric structure constant. */
+ex su3d::eval_indexed(const basic & i) const
+{
+       GINAC_ASSERT(is_of_type(i, indexed));
+       GINAC_ASSERT(i.nops() == 4);
+       GINAC_ASSERT(is_ex_of_type(i.op(0), su3d));
+
+       // Convolutions are zero
+       if (static_cast<const indexed &>(i).get_dummy_indices().size() != 0)
+               return _ex0();
+
+       // Numeric evaluation
+       if (static_cast<const indexed &>(i).all_index_values_are(info_flags::nonnegint)) {
+
+               // Sort indices
+               int v[3];
+               for (unsigned j=0; j<3; j++)
+                       v[j] = ex_to_numeric(ex_to_idx(i.op(j + 1)).get_value()).to_int();
+               if (v[0] > v[1]) std::swap(v[0], v[1]);
+               if (v[0] > v[2]) std::swap(v[0], v[2]);
+               if (v[1] > v[2]) std::swap(v[1], v[2]);
+
+#define CMPINDICES(A,B,C) ((v[0] == (A)) && (v[1] == (B)) && (v[2] == (C)))
+
+               // Check for non-zero elements
+               if (CMPINDICES(1,4,6) || CMPINDICES(1,5,7) || CMPINDICES(2,5,6)
+                || CMPINDICES(3,4,4) || CMPINDICES(3,5,5))
+                       return _ex1_2();
+               else if (CMPINDICES(2,4,7) || CMPINDICES(3,6,6) || CMPINDICES(3,7,7))
+                       return _ex_1_2();
+               else if (CMPINDICES(1,1,8) || CMPINDICES(2,2,8) || CMPINDICES(3,3,8))
+                       return sqrt(_ex3())/3;
+               else if (CMPINDICES(8,8,8))
+                       return -sqrt(_ex3())/3;
+               else if (CMPINDICES(4,4,8) || CMPINDICES(5,5,8)
+                     || CMPINDICES(6,6,8) || CMPINDICES(7,7,8))
+                       return -sqrt(_ex3())/6;
+               else
+                       return _ex0();
+       }
+
+       // No further simplifications
+       return i.hold();
+}
+
+/** Automatic symbolic evaluation of indexed antisymmetric structure constant. */
+ex su3f::eval_indexed(const basic & i) const
+{
+       GINAC_ASSERT(is_of_type(i, indexed));
+       GINAC_ASSERT(i.nops() == 4);
+       GINAC_ASSERT(is_ex_of_type(i.op(0), su3f));
+
+       // Numeric evaluation
+       if (static_cast<const indexed &>(i).all_index_values_are(info_flags::nonnegint)) {
+
+               // Sort indices, remember permutation sign
+               int v[3];
+               for (unsigned j=0; j<3; j++)
+                       v[j] = ex_to_numeric(ex_to_idx(i.op(j + 1)).get_value()).to_int();
+               int sign = 1;
+               if (v[0] > v[1]) { std::swap(v[0], v[1]); sign = -sign; }
+               if (v[0] > v[2]) { std::swap(v[0], v[2]); sign = -sign; }
+               if (v[1] > v[2]) { std::swap(v[1], v[2]); sign = -sign; }
+
+               // Check for non-zero elements
+               if (CMPINDICES(1,2,3))
+                       return sign;
+               else if (CMPINDICES(1,4,7) || CMPINDICES(2,4,6)
+                     || CMPINDICES(2,5,7) || CMPINDICES(3,4,5))
+                       return _ex1_2() * sign;
+               else if (CMPINDICES(1,5,6) || CMPINDICES(3,6,7))
+                       return _ex_1_2() * sign;
+               else if (CMPINDICES(4,5,8) || CMPINDICES(6,7,8))
+                       return sqrt(_ex3())/2 * sign;
+               else
+                       return _ex0();
+       }
+
+       // No further simplifications
+       return i.hold();
+}
+
+
+/** Contraction of an indexed symmetric structure constant with something else. */
+bool su3d::contract_with(exvector::iterator self, exvector::iterator other, exvector & v) const
+{
+       GINAC_ASSERT(is_ex_of_type(*self, indexed));
+       GINAC_ASSERT(is_ex_of_type(*other, indexed));
+       GINAC_ASSERT(self->nops() == 4);
+       GINAC_ASSERT(is_ex_of_type(self->op(0), su3d));
+
+       if (is_ex_exactly_of_type(other->op(0), su3d) || is_ex_exactly_of_type(other->op(0), su3f)) {
+
+               // Find the dummy indices of the contraction
+               exvector dummy_indices;
+               dummy_indices = ex_to_indexed(*self).get_dummy_indices(ex_to_indexed(*other));
+
+               if (is_ex_exactly_of_type(other->op(0), su3d)) {
+
+                       // d.abc*d.abc=40/3
+                       if (dummy_indices.size() == 3) {
+                               *self = numeric(40, 3);
+                               *other = _ex1();
+                               return true;
+
+                       // d.akl*d.bkl=5/3*delta.ab
+                       } else if (dummy_indices.size() == 2) {
+                               exvector a = index_set_difference(ex_to_indexed(*self).get_indices(), dummy_indices);
+                               exvector b = index_set_difference(ex_to_indexed(*other).get_indices(), dummy_indices);
+                               GINAC_ASSERT(a.size() > 0);
+                               GINAC_ASSERT(b.size() > 0);
+                               *self = numeric(5, 3) * delta_tensor(a[0], b[0]);
+                               *other = _ex1();
+                               return true;
+                       }
+
+               } else {
+
+                       // d.akl*f.bkl=0 (includes the case a=b)
+                       if (dummy_indices.size() >= 2) {
+                               *self = _ex0();
+                               *other = _ex0();
+                               return true;
+                       }
+               }
+       }
+
+       return false;
+}
+
+/** Contraction of an indexed antisymmetric structure constant with something else. */
+bool su3f::contract_with(exvector::iterator self, exvector::iterator other, exvector & v) const
+{
+       GINAC_ASSERT(is_ex_of_type(*self, indexed));
+       GINAC_ASSERT(is_ex_of_type(*other, indexed));
+       GINAC_ASSERT(self->nops() == 4);
+       GINAC_ASSERT(is_ex_of_type(self->op(0), su3f));
+
+       if (is_ex_exactly_of_type(other->op(0), su3f)) { // f*d is handled by su3d class
+
+               // Find the dummy indices of the contraction
+               exvector dummy_indices;
+               dummy_indices = ex_to_indexed(*self).get_dummy_indices(ex_to_indexed(*other));
+
+               // f.abc*f.abc=24
+               if (dummy_indices.size() == 3) {
+                       *self = 24;
+                       *other = _ex1();
+                       return true;
+
+               // f.akl*f.bkl=3*delta.ab
+               } else if (dummy_indices.size() == 2) {
+                       int sign1, sign2;
+                       ex a = permute_free_index_to_front(ex_to_indexed(*self).get_indices(), dummy_indices, sign1);
+                       ex b = permute_free_index_to_front(ex_to_indexed(*other).get_indices(), dummy_indices, sign2);
+                       *self = sign1 * sign2 * 3 * delta_tensor(a, b);
+                       *other = _ex1();
+                       return true;
+               }
+       }
+
+       return false;
+}
+

 //////////
 // global functions
 //////////
 //////////
 // global functions
 //////////


@@ -170,16 +385,31 @@ ex color_ONE(unsigned rl)
 
 ex color_T(const ex & a, unsigned rl)
 {
 
 ex color_T(const ex & a, unsigned rl)
 {

+       if (!is_ex_of_type(a, idx))
+               throw(std::invalid_argument("indices of color_T must be of type idx"));
+       if (!ex_to_idx(a).get_dim().is_equal(8))
+               throw(std::invalid_argument("index dimension for color_T must be 8"));
+

        return color(su3t(), a, rl);
 }
 
 ex color_f(const ex & a, const ex & b, const ex & c)
 {
        return color(su3t(), a, rl);
 }
 
 ex color_f(const ex & a, const ex & b, const ex & c)
 {

+       if (!is_ex_of_type(a, idx) || !is_ex_of_type(b, idx) || !is_ex_of_type(c, idx))
+               throw(std::invalid_argument("indices of color_f must be of type idx"));
+       if (!ex_to_idx(a).get_dim().is_equal(8) || !ex_to_idx(b).get_dim().is_equal(8) || !ex_to_idx(c).get_dim().is_equal(8))
+               throw(std::invalid_argument("index dimension for color_f must be 8"));
+

        return indexed(su3f(), indexed::antisymmetric, a, b, c);
 }
 
 ex color_d(const ex & a, const ex & b, const ex & c)
 {
        return indexed(su3f(), indexed::antisymmetric, a, b, c);
 }
 
 ex color_d(const ex & a, const ex & b, const ex & c)
 {

+       if (!is_ex_of_type(a, idx) || !is_ex_of_type(b, idx) || !is_ex_of_type(c, idx))
+               throw(std::invalid_argument("indices of color_d must be of type idx"));
+       if (!ex_to_idx(a).get_dim().is_equal(8) || !ex_to_idx(b).get_dim().is_equal(8) || !ex_to_idx(c).get_dim().is_equal(8))
+               throw(std::invalid_argument("index dimension for color_d must be 8"));
+

        return indexed(su3d(), indexed::symmetric, a, b, c);
 }
 
        return indexed(su3d(), indexed::symmetric, a, b, c);
 }
 

diff --git a/ginac/color.h b/ginac/color.h
index da53494644e97f95bbcb0871578184855b8a9f50..c5c57c5598942080a6f467b4ead100d17c166558 100644 (file)
--- a/ginac/color.h
+++ b/ginac/color.h
@@ -35,7 +35,8 @@ namespace GiNaC {
  *  elements from different Lie algebra representations (only objects with
  *  the same label "interact" with each other). These objects implement an
  *  abstract representation of the group, not a specific matrix
  *  elements from different Lie algebra representations (only objects with
  *  the same label "interact" with each other). These objects implement an
  *  abstract representation of the group, not a specific matrix

- *  representation. */
+ *  representation. The indices used for color objects should not have a
+ *  variance. */

 class color : public indexed
 {
        GINAC_DECLARE_REGISTERED_CLASS(color, indexed)
 class color : public indexed
 {
        GINAC_DECLARE_REGISTERED_CLASS(color, indexed)


@@ -92,8 +93,8 @@ class su3f : public tensor
        // functions overriding virtual functions from bases classes
 public:
        void print(std::ostream & os, unsigned upper_precedence=0) const;
        // functions overriding virtual functions from bases classes
 public:
        void print(std::ostream & os, unsigned upper_precedence=0) const;

-//     ex eval_indexed(const basic & i) const;
-//     bool contract_with(exvector::iterator self, exvector::iterator other, exvector & v) const;
+       ex eval_indexed(const basic & i) const;
+       bool contract_with(exvector::iterator self, exvector::iterator other, exvector & v) const;

 };
 
 /** This class represents the tensor of symmetric su(3) structure constants. */
 };
 
 /** This class represents the tensor of symmetric su(3) structure constants. */


@@ -104,8 +105,8 @@ class su3d : public tensor
        // functions overriding virtual functions from bases classes
 public:
        void print(std::ostream & os, unsigned upper_precedence=0) const;
        // functions overriding virtual functions from bases classes
 public:
        void print(std::ostream & os, unsigned upper_precedence=0) const;

-//     ex eval_indexed(const basic & i) const;
-//     bool contract_with(exvector::iterator self, exvector::iterator other, exvector & v) const;
+       ex eval_indexed(const basic & i) const;
+       bool contract_with(exvector::iterator self, exvector::iterator other, exvector & v) const;

 };
 
 
 };
 
 

diff --git a/ginac/idx.cpp b/ginac/idx.cpp
index 994d532ecc80c1da10729b58ffadeefc87f35fee..d8e8e70c6368b4a3f6116f1cfbe0d58be1f6a384 100644 (file)
--- a/ginac/idx.cpp
+++ b/ginac/idx.cpp
@@ -305,4 +305,90 @@ bool is_dummy_pair(const ex & e1, const ex & e2)
        return is_dummy_pair(ex_to_idx(e1), ex_to_idx(e2));
 }
 
        return is_dummy_pair(ex_to_idx(e1), ex_to_idx(e2));
 }
 

+/** Bring a vector of indices into a canonic order. Dummy indices will lie
+ *  next to each other after the sorting. */
+static void sort_index_vector(exvector &v)
+{
+       // Nothing to sort if less than 2 elements
+       if (v.size() < 2)
+               return;
+
+       // Simple bubble sort algorithm should be sufficient for the small
+       // number of indices expected
+       exvector::iterator it1 = v.begin(), itend = v.end(), next_to_last_idx = itend - 1;
+       while (it1 != next_to_last_idx) {
+               exvector::iterator it2 = it1 + 1;
+               while (it2 != itend) {
+                       if (it1->compare(*it2) > 0)
+                               it1->swap(*it2);
+                       it2++;
+               }
+               it1++;
+       }
+}
+
+
+void find_free_and_dummy(exvector::const_iterator it, exvector::const_iterator itend, exvector & out_free, exvector & out_dummy)
+{
+       out_free.clear();
+       out_dummy.clear();
+
+       // No indices? Then do nothing
+       if (it == itend)
+               return;
+
+       // Only one index? Then it is a free one if it's not numeric
+       if (itend - it == 1) {
+               if (ex_to_idx(*it).is_symbolic())
+                       out_free.push_back(*it);
+               return;
+       }
+
+       // Sort index vector. This will cause dummy indices come to lie next
+       // to each other (because the sort order is defined to guarantee this).
+       exvector v(it, itend);
+       sort_index_vector(v);
+
+       // Find dummy pairs and free indices
+       it = v.begin(); itend = v.end();
+       exvector::const_iterator last = it++;
+       while (it != itend) {
+               if (is_dummy_pair(*it, *last)) {
+                       out_dummy.push_back(*last);
+                       it++;
+                       if (it == itend)
+                               return;
+               } else {
+                       if (!it->is_equal(*last) && ex_to_idx(*last).is_symbolic())
+                               out_free.push_back(*last);
+               }
+               last = it++;
+       }
+       if (ex_to_idx(*last).is_symbolic())
+               out_free.push_back(*last);
+}
+
+exvector index_set_difference(const exvector & set1, const exvector & set2)
+{
+       exvector ret;
+
+       exvector::const_iterator ait = set1.begin(), aitend = set1.end();
+       while (ait != aitend) {
+               exvector::const_iterator bit = set2.begin(), bitend = set2.end();
+               bool found = false;
+               while (bit != bitend) {
+                       if (ait->is_equal(*bit)) {
+                               found = true;
+                               break;
+                       }
+                       bit++;
+               }
+               if (!found)
+                       ret.push_back(*ait);
+               ait++;
+       }
+
+       return ret;
+}
+

 } // namespace GiNaC
 } // namespace GiNaC

diff --git a/ginac/idx.h b/ginac/idx.h
index ea885feef1d56c86f043b328b5b4f279c010fbe7..aa7f94e1fdf3fa89ae8846e96de0b04db5855854 100644 (file)
--- a/ginac/idx.h
+++ b/ginac/idx.h
@@ -141,6 +141,57 @@ bool is_dummy_pair(const idx & i1, const idx & i2);
 /** Check whether two expressions form a dummy index pair. */
 bool is_dummy_pair(const ex & e1, const ex & e2);
 
 /** Check whether two expressions form a dummy index pair. */
 bool is_dummy_pair(const ex & e1, const ex & e2);
 

+/** Given a vector of indices, split them into two vectors, one containing
+ *  the free indices, the other containing the dummy indices (numeric
+ *  indices are neither free nor dummy ones).
+ *
+ *  @param it Pointer to start of index vector
+ *  @param itend Pointer to end of index vector
+ *  @param out_free Vector of free indices (returned, sorted)
+ *  @param out_dummy Vector of dummy indices (returned, sorted) */
+void find_free_and_dummy(exvector::const_iterator it, exvector::const_iterator itend, exvector & out_free, exvector & out_dummy);
+
+/** Given a vector of indices, split them into two vectors, one containing
+ *  the free indices, the other containing the dummy indices (numeric
+ *  indices are neither free nor dummy ones).
+ *
+ *  @param v Index vector
+ *  @param out_free Vector of free indices (returned, sorted)
+ *  @param out_dummy Vector of dummy indices (returned, sorted) */
+inline void find_free_and_dummy(const exvector & v, exvector & out_free, exvector & out_dummy)
+{
+       find_free_and_dummy(v.begin(), v.end(), out_free, out_dummy);
+}
+
+/** Given a vector of indices, find the dummy indices.
+ *
+ *  @param v Index vector
+ *  @param out_dummy Vector of dummy indices (returned, sorted) */
+inline void find_dummy_indices(const exvector & v, exvector & out_dummy)
+{
+       exvector free_indices;
+       find_free_and_dummy(v.begin(), v.end(), free_indices, out_dummy);
+}
+
+/** Count the number of dummy index pairs in an index vector. */
+inline unsigned count_dummy_indices(const exvector & v)
+{
+       exvector free_indices, dummy_indices;
+       find_free_and_dummy(v.begin(), v.end(), free_indices, dummy_indices);
+       return dummy_indices.size();
+}
+
+/** Count the number of dummy index pairs in an index vector. */
+inline unsigned count_free_indices(const exvector & v)
+{
+       exvector free_indices, dummy_indices;
+       find_free_and_dummy(v.begin(), v.end(), free_indices, dummy_indices);
+       return free_indices.size();
+}
+
+/** Given two index vectors, find those indices that appear in the first
+ *  vector but not in the second one (asymmetric set difference). */
+exvector index_set_difference(const exvector & set1, const exvector & set2);

 
 } // namespace GiNaC
 
 
 } // namespace GiNaC
 

diff --git a/ginac/indexed.cpp b/ginac/indexed.cpp
index dfd9a3812876c9707e9e35a0f6ae03e1dddddf15..0ab989e801baccc8a94d60a7df7306da1822a78f 100644 (file)
--- a/ginac/indexed.cpp
+++ b/ginac/indexed.cpp
@@ -445,48 +445,6 @@ void indexed::assert_all_indices_of_type_idx(void) const
 // global functions
 //////////
 
 // global functions
 //////////
 

-/** Given a vector of indices, split them into two vectors, one containing
- *  the free indices, the other containing the dummy indices. */
-static void find_free_and_dummy(exvector::const_iterator it, exvector::const_iterator itend, exvector & out_free, exvector & out_dummy)
-{
-       out_free.clear();
-       out_dummy.clear();
-
-       // No indices? Then do nothing
-       if (it == itend)
-               return;
-
-       // Only one index? Then it is a free one if it's not numeric
-       if (itend - it == 1) {
-               if (ex_to_idx(*it).is_symbolic())
-                       out_free.push_back(*it);
-               return;
-       }
-
-       // Sort index vector. This will cause dummy indices come to lie next
-       // to each other (because the sort order is defined to guarantee this).
-       exvector v(it, itend);
-       sort_index_vector(v);
-
-       // Find dummy pairs and free indices
-       it = v.begin(); itend = v.end();
-       exvector::const_iterator last = it++;
-       while (it != itend) {
-               if (is_dummy_pair(*it, *last)) {
-                       out_dummy.push_back(*last);
-                       it++;
-                       if (it == itend)
-                               return;
-               } else {
-                       if (!it->is_equal(*last) && ex_to_idx(*last).is_symbolic())
-                               out_free.push_back(*last);
-               }
-               last = it++;
-       }
-       if (ex_to_idx(*last).is_symbolic())
-               out_free.push_back(*last);
-}
-

 /** Check whether two sorted index vectors are consistent (i.e. equal). */
 static bool indices_consistent(const exvector & v1, const exvector & v2)
 {
 /** Check whether two sorted index vectors are consistent (i.e. equal). */
 static bool indices_consistent(const exvector & v1, const exvector & v2)
 {


@@ -505,6 +463,12 @@ static bool indices_consistent(const exvector & v1, const exvector & v2)
        return true;
 }
 
        return true;
 }
 

+exvector indexed::get_indices(void) const
+{
+       GINAC_ASSERT(seq.size() >= 1);
+       return exvector(seq.begin() + 1, seq.end());
+}
+

 exvector indexed::get_dummy_indices(void) const
 {
        exvector free_indices, dummy_indices;
 exvector indexed::get_dummy_indices(void) const
 {
        exvector free_indices, dummy_indices;


@@ -512,6 +476,16 @@ exvector indexed::get_dummy_indices(void) const
        return dummy_indices;
 }
 
        return dummy_indices;
 }
 

+exvector indexed::get_dummy_indices(const indexed & other) const
+{
+       exvector indices = get_free_indices();
+       exvector other_indices = other.get_free_indices();
+       indices.insert(indices.end(), other_indices.begin(), other_indices.end());
+       exvector dummy_indices;
+       find_dummy_indices(indices, dummy_indices);
+       return dummy_indices;
+}
+

 exvector indexed::get_free_indices(void) const
 {
        exvector free_indices, dummy_indices;
 exvector indexed::get_free_indices(void) const
 {
        exvector free_indices, dummy_indices;


@@ -545,7 +519,7 @@ exvector mul::get_free_indices(void) const
 
        // And remove the dummy indices
        exvector free_indices, dummy_indices;
 
        // And remove the dummy indices
        exvector free_indices, dummy_indices;

-       find_free_and_dummy(un.begin(), un.end(), free_indices, dummy_indices);
+       find_free_and_dummy(un, free_indices, dummy_indices);

        return free_indices;
 }
 
        return free_indices;
 }
 


@@ -560,7 +534,7 @@ exvector ncmul::get_free_indices(void) const
 
        // And remove the dummy indices
        exvector free_indices, dummy_indices;
 
        // And remove the dummy indices
        exvector free_indices, dummy_indices;

-       find_free_and_dummy(un.begin(), un.end(), free_indices, dummy_indices);
+       find_free_and_dummy(un, free_indices, dummy_indices);

        return free_indices;
 }
 
        return free_indices;
 }
 


@@ -596,7 +570,7 @@ ex simplify_indexed_product(const ex & e, exvector & free_indices, const scalar_
                        } else if (is_ex_exactly_of_type(f, ncmul)) {
                                // Noncommutative factor found, split it as well
                                non_commutative = true; // everything becomes noncommutative, ncmul will sort out the commutative factors later
                        } else if (is_ex_exactly_of_type(f, ncmul)) {
                                // Noncommutative factor found, split it as well
                                non_commutative = true; // everything becomes noncommutative, ncmul will sort out the commutative factors later

-                               for (int j=0; j<f.nops(); i++)
+                               for (int j=0; j<f.nops(); j++)

                                        v.push_back(f.op(j));
                        } else
                                v.push_back(f);
                                        v.push_back(f.op(j));
                        } else
                                v.push_back(f);


@@ -624,7 +598,7 @@ try_again:
                        exvector un(ex_to_indexed(*it1).seq.begin() + 1, ex_to_indexed(*it1).seq.end());
                        un.insert(un.end(), ex_to_indexed(*it2).seq.begin() + 1, ex_to_indexed(*it2).seq.end());
                        exvector free, dummy;
                        exvector un(ex_to_indexed(*it1).seq.begin() + 1, ex_to_indexed(*it1).seq.end());
                        un.insert(un.end(), ex_to_indexed(*it2).seq.begin() + 1, ex_to_indexed(*it2).seq.end());
                        exvector free, dummy;

-                       find_free_and_dummy(un.begin(), un.end(), free, dummy);
+                       find_free_and_dummy(un, free, dummy);

                        if (dummy.size() == 0)
                                continue;
 
                        if (dummy.size() == 0)
                                continue;
 


@@ -677,7 +651,7 @@ try_again:
                }
                it1++;
        }
                }
                it1++;
        }

-       find_free_and_dummy(un.begin(), un.end(), free_indices, dummy_indices);
+       find_free_and_dummy(un, free_indices, dummy_indices);

 
        ex r;
        if (something_changed) {
 
        ex r;
        if (something_changed) {

diff --git a/ginac/indexed.h b/ginac/indexed.h
index b06fad1040a702a8354248f3ef68834297c88cc8..301f2d748669091813efc21b3011f8ca7110ac2c 100644 (file)
--- a/ginac/indexed.h
+++ b/ginac/indexed.h
@@ -173,9 +173,16 @@ public:
         *  @see class info_flags */
        bool all_index_values_are(unsigned inf) const;
 
         *  @see class info_flags */
        bool all_index_values_are(unsigned inf) const;
 

+       /** Return a vector containing the object's indices. */
+       exvector get_indices(void) const;
+

        /** Return a vector containing the dummy indices of the object, if any. */
        exvector get_dummy_indices(void) const;
 
        /** Return a vector containing the dummy indices of the object, if any. */
        exvector get_dummy_indices(void) const;
 

+       /** Return a vector containing the dummy indices in the contraction with
+        *  another indexed object. */
+       exvector get_dummy_indices(const indexed & other) const;
+

 protected:
        void printrawindices(std::ostream & os) const;
        void printtreeindices(std::ostream & os, unsigned indent) const;
 protected:
        void printrawindices(std::ostream & os) const;
        void printtreeindices(std::ostream & os, unsigned indent) const;

diff --git a/ginac/tensor.cpp b/ginac/tensor.cpp
index ed6ad9dc4fbe373ac590d40b6ab48e64d424cbf2..b180ec9e9a7e6966cdabb06ddbdd686e201aa479 100644 (file)
--- a/ginac/tensor.cpp
+++ b/ginac/tensor.cpp
@@ -179,6 +179,15 @@ ex tensdelta::eval_indexed(const basic & i) const
        if (is_dummy_pair(i1, i2))
                return i1.get_dim();
 
        if (is_dummy_pair(i1, i2))
                return i1.get_dim();
 

+       // Numeric evaluation
+       if (static_cast<const indexed &>(i).all_index_values_are(info_flags::integer)) {
+               int n1 = ex_to_numeric(i1.get_value()).to_int(), n2 = ex_to_numeric(i2.get_value()).to_int();
+               if (n1 == n2)
+                       return _ex1();
+               else
+                       return _ex0();
+       }
+

        // No further simplifications
        return i.hold();
 }
        // No further simplifications
        return i.hold();
 }

diff --git a/ginac/tinfos.h b/ginac/tinfos.h
index 767285576bfb19eba4ce71ffc291621e7eb910c2..d71f251d09ca1a9f2b324d2705c79dce70d7f6b3 100644 (file)
--- a/ginac/tinfos.h
+++ b/ginac/tinfos.h
@@ -73,6 +73,7 @@ const unsigned TINFO_su3t          = 0x000e1009U;
 const unsigned TINFO_su3f          = 0x000e100aU;
 const unsigned TINFO_su3d          = 0x000e100bU;
 const unsigned TINFO_diracgamma    = 0x000e100cU;
 const unsigned TINFO_su3f          = 0x000e100aU;
 const unsigned TINFO_su3d          = 0x000e100bU;
 const unsigned TINFO_diracgamma    = 0x000e100cU;

+const unsigned TINFO_diracone      = 0x000e100dU;

 
 } // namespace GiNaC
 
 
 } // namespace GiNaC