- simplify_indexed() renames dummy indices so, e.g., "a.i*a.i+a.j*a.j" gets
authorChristian Bauer <Christian.Bauer@uni-mainz.de>
Fri, 18 May 2001 00:53:07 +0000 (00:53 +0000)
committerChristian Bauer <Christian.Bauer@uni-mainz.de>
Fri, 18 May 2001 00:53:07 +0000 (00:53 +0000)
  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()

ginac/clifford.cpp
ginac/clifford.h
ginac/expairseq.cpp
ginac/indexed.cpp
ginac/indexed.h
ginac/ncmul.cpp
ginac/ncmul.h
ginac/normal.cpp

index 5c7f7c3..a903881 100644 (file)
@@ -446,4 +446,62 @@ ex dirac_trace(const ex & e, unsigned char rl, const ex & trONE)
        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
index 42be538..487c12d 100644 (file)
@@ -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);
 
+/** 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
 
index 9f57061..b6a0f37 100644 (file)
@@ -991,7 +991,6 @@ void expairseq::make_flat(const epvector &v)
        return;
 }
 
-
 /** Brings this expairseq into a sorted (canonical) form. */
 void expairseq::canonicalize(void)
 {
index 3240c8d..435d952 100644 (file)
@@ -21,6 +21,7 @@
  */
 
 #include <stdexcept>
+#include <algorithm>
 
 #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.
 
-/** 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.
@@ -570,9 +546,75 @@ exvector power::get_free_indices(void) const
        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. */
-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)
@@ -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)));
-                                       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
@@ -690,14 +732,14 @@ contraction_done:
        }
 
        // 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++;
        }
-       find_free_and_dummy(un, free_indices, dummy_indices);
+       find_free_and_dummy(un, free_indices, local_dummy_indices);
 
        ex r;
        if (something_changed)
@@ -705,6 +747,9 @@ contraction_done:
        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))
@@ -714,17 +759,18 @@ contraction_done:
 }
 
 /** 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
+       // (and perform dummy index renaming if 
        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
@@ -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;
-                       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;
@@ -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())))
-               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();
@@ -769,15 +815,15 @@ ex simplify_indexed(const ex & e, exvector & free_indices, const scalar_products
 
 ex simplify_indexed(const ex & e)
 {
-       exvector free_indices;
+       exvector free_indices, dummy_indices;
        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)
 {
-       exvector free_indices;
-       return simplify_indexed(e, free_indices, sp);
+       exvector free_indices, dummy_indices;
+       return simplify_indexed(e, free_indices, dummy_indices, sp);
 }
 
 //////////
index 0147185..26bf9a5 100644 (file)
@@ -39,8 +39,8 @@ class indexed : public 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:
index e775d82..abe72ef 100644 (file)
@@ -432,11 +432,6 @@ ex ncmul::eval(int level) const
                                                                                  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);
index 05f9994..05b96da 100644 (file)
@@ -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 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;
index 369e2c6..d3a407d 100644 (file)
@@ -36,7 +36,6 @@
 #include "inifcns.h"
 #include "lst.h"
 #include "mul.h"
-#include "ncmul.h"
 #include "numeric.h"
 #include "power.h"
 #include "relational.h"