- simplify_indexed() renames dummy indices so, e.g., "a.i*a.i+a.j*a.j" gets

  simplified to "2*a.i*a.i" (or 2*a.j*a.j, you can't know which one)
- fixed possible crash when calling subs() on expressions with non-commutative
  products
- added canonicalize_clifford()

diff --git a/ginac/clifford.cpp b/ginac/clifford.cpp
index 5c7f7c3b1623e2941131680ed8aa8463f9de2532..a9038818d3e618c28457f3a2b75a6131578a911c 100644 (file)
--- a/ginac/clifford.cpp
+++ b/ginac/clifford.cpp
@@ -446,4 +446,62 @@ ex dirac_trace(const ex & e, unsigned char rl, const ex & trONE)
        return _ex0();
 }
 
        return _ex0();
 }
 

+ex canonicalize_clifford(const ex & e)
+{
+       if (is_ex_exactly_of_type(e, add)) {
+
+               ex sum = _ex0();
+               for (unsigned i=0; i<e.nops(); i++)
+                       sum += canonicalize_clifford(e.op(i));
+               return sum;
+
+       } else if (is_ex_exactly_of_type(e, mul)) {
+
+               ex prod = _ex1();
+               for (unsigned i=0; i<e.nops(); i++)
+                       prod *= canonicalize_clifford(e.op(i));
+               return prod;
+
+       } else if (is_ex_exactly_of_type(e, ncmul)) {
+
+               // Expand product, if necessary
+               ex e_expanded = e.expand();
+               if (!is_ex_of_type(e_expanded, ncmul))
+                       return canonicalize_clifford(e_expanded);
+
+               if (!is_ex_of_type(e.op(0), clifford))
+                       return e;
+
+               exvector v;
+               v.reserve(e.nops());
+               for (int i=0; i<e.nops(); i++)
+                       v.push_back(e.op(i));
+
+               // Stupid bubble sort because we only want to swap adjacent gammas
+               exvector::iterator itstart = v.begin(), itend = v.end(), next_to_last = itend - 1;
+               if (is_ex_of_type(itstart->op(0), diracgamma5))
+                       itstart++;
+               while (next_to_last != itstart) {
+                       exvector::iterator it = itstart;
+                       while (it != next_to_last) {
+                               if (it[0].op(1).compare(it[1].op(1)) > 0) {
+                                       ex save0 = it[0], save1 = it[1];
+                                       it[0] = lorentz_g(it[0].op(1), it[1].op(1));
+                                       it[1] = _ex2();
+                                       ex sum = ncmul(v);
+                                       it[0] = save1;
+                                       it[1] = save0;
+                                       sum -= ncmul(v);
+                                       return canonicalize_clifford(sum);
+                               }
+                               it++;
+                       }
+                       next_to_last--;
+               }
+               return ncmul(v);
+       }
+
+       return e;
+}
+

 } // namespace GiNaC
 } // namespace GiNaC

diff --git a/ginac/clifford.h b/ginac/clifford.h
index 42be538f5f5bdbe9beb29dd95d67af6ab951b792..487c12d7f7d1e539d8252ea6a7ffcc98093bfeea 100644 (file)
--- a/ginac/clifford.h
+++ b/ginac/clifford.h
@@ -149,6 +149,10 @@ ex dirac_slash(const ex & e, const ex & dim, unsigned char rl = 0);
  *  @param trONE Expression to be returned as the trace of the unit matrix */
 ex dirac_trace(const ex & e, unsigned char rl = 0, const ex & trONE = 4);
 
  *  @param trONE Expression to be returned as the trace of the unit matrix */
 ex dirac_trace(const ex & e, unsigned char rl = 0, const ex & trONE = 4);
 

+/** Bring all products of clifford objects in an expression into a canonical
+ *  order. This is not necessarily the most simple form but it will allow
+ *  to checking two expressions for equality. */
+ex canonicalize_clifford(const ex & e);

 
 } // namespace GiNaC
 
 
 } // namespace GiNaC
 

diff --git a/ginac/expairseq.cpp b/ginac/expairseq.cpp
index 9f570617ac7b6185edefecc003142cfa7a524e48..b6a0f372adb8e762abb4d7ff541b53b54b6d6ddb 100644 (file)
--- a/ginac/expairseq.cpp
+++ b/ginac/expairseq.cpp
@@ -991,7 +991,6 @@ void expairseq::make_flat(const epvector &v)
        return;
 }
 
        return;
 }
 

-

 /** Brings this expairseq into a sorted (canonical) form. */
 void expairseq::canonicalize(void)
 {
 /** Brings this expairseq into a sorted (canonical) form. */
 void expairseq::canonicalize(void)
 {

diff --git a/ginac/indexed.cpp b/ginac/indexed.cpp
index 3240c8dd7fe86293918c0aabf807cecd91e101b8..435d952784707416c729848ec6914a46fdb43b72 100644 (file)
--- a/ginac/indexed.cpp
+++ b/ginac/indexed.cpp
@@ -21,6 +21,7 @@
  */
 
 #include <stdexcept>
  */
 
 #include <stdexcept>

+#include <algorithm>

 
 #include "indexed.h"
 #include "idx.h"
 
 #include "indexed.h"
 #include "idx.h"


@@ -253,31 +254,6 @@ int indexed::compare_same_type(const basic & other) const
 // reorder index pairs with known symmetry properties, while sort_index_vector()
 // always sorts the whole vector.
 
 // reorder index pairs with known symmetry properties, while sort_index_vector()
 // always sorts the whole vector.
 

-/** Bring a vector of indices into a canonic order (don't care about the
- *  symmetry of the objects carrying the indices). Dummy indices will lie
- *  next to each other after the sorting.
- *
- *  @param v Index vector to be sorted */
-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++;
-       }
-}
-

 /** Bring a vector of indices into a canonic order. This operation only makes
  *  sense if the object carrying these indices is either symmetric or totally
  *  antisymmetric with respect to the indices.
 /** Bring a vector of indices into a canonic order. This operation only makes
  *  sense if the object carrying these indices is either symmetric or totally
  *  antisymmetric with respect to the indices.


@@ -570,9 +546,75 @@ exvector power::get_free_indices(void) const
        return basis.get_free_indices();
 }
 
        return basis.get_free_indices();
 }
 

+/* Function object for STL sort() */
+struct ex_is_less {
+       bool operator() (const ex &lh, const ex &rh) const
+       {
+               return lh.compare(rh) < 0;
+       }
+};
+
+/** Rename dummy indices in an expression.
+ *
+ *  @param e Expression to be worked on
+ *  @param local_dummy_indices The set of dummy indices that appear in the
+ *    expression "e"
+ *  @param global_dummy_indices The set of dummy indices that have appeared
+ *    before and which we would like to use in "e", too. This gets updated
+ *    by the function */
+static ex rename_dummy_indices(const ex & e, exvector & global_dummy_indices, exvector & local_dummy_indices)
+{
+       int global_size = global_dummy_indices.size(),
+           local_size = local_dummy_indices.size();
+
+       // Any local dummy indices at all?
+       if (local_size == 0)
+               return e;
+
+       sort(local_dummy_indices.begin(), local_dummy_indices.end(), ex_is_less());
+
+       if (global_size < local_size) {
+
+               // More local indices than we encountered before, add the new ones
+               // to the global set
+               int remaining = local_size - global_size;
+               exvector::const_iterator it = local_dummy_indices.begin(), itend = local_dummy_indices.end();
+               while (it != itend && remaining > 0) {
+                       exvector::const_iterator git = global_dummy_indices.begin(), gitend = global_dummy_indices.end();
+                       while (git != gitend) {
+                               if (it->is_equal(*git))
+                                       goto found;
+                               git++;
+                       }
+                       global_dummy_indices.push_back(*it);
+                       global_size++;
+                       remaining--;
+found:         it++;
+               }
+               sort(global_dummy_indices.begin(), global_dummy_indices.end(), ex_is_less());
+       }
+
+       // Replace index symbols in expression
+       GINAC_ASSERT(local_size <= global_size);
+       bool all_equal = true;
+       lst local_syms, global_syms;
+       for (unsigned i=0; i<local_size; i++) {
+               ex loc_sym = local_dummy_indices[i].op(0);
+               ex glob_sym = global_dummy_indices[i].op(0);
+               if (!loc_sym.is_equal(glob_sym))
+                       all_equal = false;
+               local_syms.append(loc_sym);
+               global_syms.append(glob_sym);
+       }
+       if (all_equal)
+               return e;
+       else
+               return e.subs(local_syms, global_syms);
+}
+

 /** Simplify product of indexed expressions (commutative, noncommutative and
  *  simple squares), return list of free indices. */
 /** Simplify product of indexed expressions (commutative, noncommutative and
  *  simple squares), return list of free indices. */

-ex simplify_indexed_product(const ex & e, exvector & free_indices, const scalar_products & sp)
+ex simplify_indexed_product(const ex & e, exvector & free_indices, exvector & dummy_indices, const scalar_products & sp)

 {
        // Remember whether the product was commutative or noncommutative
        // (because we chop it into factors and need to reassemble later)
 {
        // Remember whether the product was commutative or noncommutative
        // (because we chop it into factors and need to reassemble later)


@@ -677,7 +719,7 @@ contraction_done:
                                        // simplify_ncmul() the chance to re-order and canonicalize
                                        // the product
                                        ex r = (non_commutative ? ex(ncmul(v)) : ex(mul(v)));
                                        // simplify_ncmul() the chance to re-order and canonicalize
                                        // the product
                                        ex r = (non_commutative ? ex(ncmul(v)) : ex(mul(v)));

-                                       return simplify_indexed(r, free_indices, sp);
+                                       return simplify_indexed(r, free_indices, dummy_indices, sp);

                                }
 
                                // Both objects may have new indices now or they might
                                }
 
                                // Both objects may have new indices now or they might


@@ -690,14 +732,14 @@ contraction_done:
        }
 
        // Find free indices (concatenate them all and call find_free_and_dummy())
        }
 
        // Find free indices (concatenate them all and call find_free_and_dummy())

-       exvector un, dummy_indices;
+       exvector un, local_dummy_indices;

        it1 = v.begin(); itend = v.end();
        while (it1 != itend) {
                exvector free_indices_of_factor = it1->get_free_indices();
                un.insert(un.end(), free_indices_of_factor.begin(), free_indices_of_factor.end());
                it1++;
        }
        it1 = v.begin(); itend = v.end();
        while (it1 != itend) {
                exvector free_indices_of_factor = it1->get_free_indices();
                un.insert(un.end(), free_indices_of_factor.begin(), free_indices_of_factor.end());
                it1++;
        }

-       find_free_and_dummy(un, free_indices, dummy_indices);
+       find_free_and_dummy(un, free_indices, local_dummy_indices);

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


@@ -705,6 +747,9 @@ contraction_done:
        else
                r = e;
 
        else
                r = e;
 

+       // Dummy index renaming
+       r = rename_dummy_indices(r, dummy_indices, local_dummy_indices);
+

        // Product of indexed object with a scalar?
        if (is_ex_exactly_of_type(r, mul) && r.nops() == 2
         && is_ex_exactly_of_type(r.op(1), numeric) && is_ex_of_type(r.op(0), indexed))
        // Product of indexed object with a scalar?
        if (is_ex_exactly_of_type(r, mul) && r.nops() == 2
         && is_ex_exactly_of_type(r.op(1), numeric) && is_ex_of_type(r.op(0), indexed))


@@ -714,17 +759,18 @@ contraction_done:
 }
 
 /** Simplify indexed expression, return list of free indices. */
 }
 
 /** Simplify indexed expression, return list of free indices. */

-ex simplify_indexed(const ex & e, exvector & free_indices, const scalar_products & sp)
+ex simplify_indexed(const ex & e, exvector & free_indices, exvector & dummy_indices, const scalar_products & sp)

 {
        // Expand the expression
        ex e_expanded = e.expand();
 
        // Simplification of single indexed object: just find the free indices
 {
        // Expand the expression
        ex e_expanded = e.expand();
 
        // Simplification of single indexed object: just find the free indices

+       // (and perform dummy index renaming if 

        if (is_ex_of_type(e_expanded, indexed)) {
                const indexed &i = ex_to_indexed(e_expanded);
        if (is_ex_of_type(e_expanded, indexed)) {
                const indexed &i = ex_to_indexed(e_expanded);

-               exvector dummy_indices;
-               find_free_and_dummy(i.seq.begin() + 1, i.seq.end(), free_indices, dummy_indices);
-               return e_expanded;
+               exvector local_dummy_indices;
+               find_free_and_dummy(i.seq.begin() + 1, i.seq.end(), free_indices, local_dummy_indices);
+               return rename_dummy_indices(e_expanded, dummy_indices, local_dummy_indices);

        }
 
        // Simplification of sum = sum of simplifications, check consistency of
        }
 
        // Simplification of sum = sum of simplifications, check consistency of


@@ -736,7 +782,7 @@ ex simplify_indexed(const ex & e, exvector & free_indices, const scalar_products
 
                for (unsigned i=0; i<e_expanded.nops(); i++) {
                        exvector free_indices_of_term;
 
                for (unsigned i=0; i<e_expanded.nops(); i++) {
                        exvector free_indices_of_term;

-                       ex term = simplify_indexed(e_expanded.op(i), free_indices_of_term, sp);
+                       ex term = simplify_indexed(e_expanded.op(i), free_indices_of_term, dummy_indices, sp);

                        if (!term.is_zero()) {
                                if (first) {
                                        free_indices = free_indices_of_term;
                        if (!term.is_zero()) {
                                if (first) {
                                        free_indices = free_indices_of_term;


@@ -760,7 +806,7 @@ ex simplify_indexed(const ex & e, exvector & free_indices, const scalar_products
        if (is_ex_exactly_of_type(e_expanded, mul)
         || is_ex_exactly_of_type(e_expanded, ncmul)
         || (is_ex_exactly_of_type(e_expanded, power) && is_ex_of_type(e_expanded.op(0), indexed) && e_expanded.op(1).is_equal(_ex2())))
        if (is_ex_exactly_of_type(e_expanded, mul)
         || is_ex_exactly_of_type(e_expanded, ncmul)
         || (is_ex_exactly_of_type(e_expanded, power) && is_ex_of_type(e_expanded.op(0), indexed) && e_expanded.op(1).is_equal(_ex2())))

-               return simplify_indexed_product(e_expanded, free_indices, sp);
+               return simplify_indexed_product(e_expanded, free_indices, dummy_indices, sp);

 
        // Cannot do anything
        free_indices.clear();
 
        // Cannot do anything
        free_indices.clear();


@@ -769,15 +815,15 @@ ex simplify_indexed(const ex & e, exvector & free_indices, const scalar_products
 
 ex simplify_indexed(const ex & e)
 {
 
 ex simplify_indexed(const ex & e)
 {

-       exvector free_indices;
+       exvector free_indices, dummy_indices;

        scalar_products sp;
        scalar_products sp;

-       return simplify_indexed(e, free_indices, sp);
+       return simplify_indexed(e, free_indices, dummy_indices, sp);

 }
 
 ex simplify_indexed(const ex & e, const scalar_products & sp)
 {
 }
 
 ex simplify_indexed(const ex & e, const scalar_products & sp)
 {

-       exvector free_indices;
-       return simplify_indexed(e, free_indices, sp);
+       exvector free_indices, dummy_indices;
+       return simplify_indexed(e, free_indices, dummy_indices, sp);

 }
 
 //////////
 }
 
 //////////

diff --git a/ginac/indexed.h b/ginac/indexed.h
index 0147185ddd4630183ee7e78b23db134015fccadd..26bf9a5eec1c2563410ea9195b38b9ef1bb8beff 100644 (file)
--- a/ginac/indexed.h
+++ b/ginac/indexed.h
@@ -39,8 +39,8 @@ class indexed : public exprseq
 {
        GINAC_DECLARE_REGISTERED_CLASS(indexed, exprseq)
 
 {
        GINAC_DECLARE_REGISTERED_CLASS(indexed, exprseq)
 

-       friend ex simplify_indexed(const ex & e, exvector & free_indices, const scalar_products & sp);
-       friend ex simplify_indexed_product(const ex & e, exvector & free_indices, const scalar_products & sp);
+       friend ex simplify_indexed(const ex & e, exvector & free_indices, exvector & dummy_indices, const scalar_products & sp);
+       friend ex simplify_indexed_product(const ex & e, exvector & free_indices, exvector & dummy_indices, const scalar_products & sp);

 
        // types
 public:
 
        // types
 public:

diff --git a/ginac/ncmul.cpp b/ginac/ncmul.cpp
index e775d82e9e2d85207dd1330bd494d801d195ab9e..abe72ef1f726e64e560ebaca45061d3f4ddc04e9 100644 (file)
--- a/ginac/ncmul.cpp
+++ b/ginac/ncmul.cpp
@@ -432,11 +432,6 @@ ex ncmul::eval(int level) const
                                                                                  status_flags::evaluated);
 }
 
                                                                                  status_flags::evaluated);
 }
 

-ex ncmul::subs(const lst & ls, const lst & lr) const
-{
-       return ncmul(subschildren(ls, lr));
-}
-

 ex ncmul::thisexprseq(const exvector & v) const
 {
        return (new ncmul(v))->setflag(status_flags::dynallocated);
 ex ncmul::thisexprseq(const exvector & v) const
 {
        return (new ncmul(v))->setflag(status_flags::dynallocated);

diff --git a/ginac/ncmul.h b/ginac/ncmul.h
index 05f9994929ad0fd28263212c0e224b129fec1e51..05b96da375d2e72936195ba9bef5976367d1307d 100644 (file)
--- a/ginac/ncmul.h
+++ b/ginac/ncmul.h
@@ -60,7 +60,6 @@ public:
        ex expand(unsigned options=0) const;
        ex coeff(const ex & s, int n=1) const;
        ex eval(int level=0) const;
        ex expand(unsigned options=0) const;
        ex coeff(const ex & s, int n=1) const;
        ex eval(int level=0) const;

-       ex subs(const lst & ls, const lst & lr) const;

        exvector get_free_indices(void) const;
        ex thisexprseq(const exvector & v) const;
        ex thisexprseq(exvector * vp) const;
        exvector get_free_indices(void) const;
        ex thisexprseq(const exvector & v) const;
        ex thisexprseq(exvector * vp) const;

diff --git a/ginac/normal.cpp b/ginac/normal.cpp
index 369e2c656449ff296385397b1bce19c1b1b356af..d3a407d50e713ba633054784f263976360f5dc49 100644 (file)
--- a/ginac/normal.cpp
+++ b/ginac/normal.cpp
@@ -36,7 +36,6 @@
 #include "inifcns.h"
 #include "lst.h"
 #include "mul.h"
 #include "inifcns.h"
 #include "lst.h"
 #include "mul.h"

-#include "ncmul.h"

 #include "numeric.h"
 #include "power.h"
 #include "relational.h"
 #include "numeric.h"
 #include "power.h"
 #include "relational.h"