- dummy index renamer didn't account for internal dummy indices of objects
authorChristian Bauer <Christian.Bauer@uni-mainz.de>
Sat, 19 May 2001 00:42:07 +0000 (00:42 +0000)
committerChristian Bauer <Christian.Bauer@uni-mainz.de>
Sat, 19 May 2001 00:42:07 +0000 (00:42 +0000)
  in products; a~mu.mu-a~nu.nu gets simplified to 0 now
- made a little more use of STL facilities for exvectors, especially in the
  indexed stuff; append_exvector_to_exvector() and index_set_difference() are
  gone and utils.h defines the functors ex_is_less and ex_is_equal

ginac/basic.h
ginac/clifford.cpp
ginac/color.cpp
ginac/idx.cpp
ginac/idx.h
ginac/indexed.cpp
ginac/matrix.h
ginac/power.cpp
ginac/registrar.h
ginac/utils.cpp
ginac/utils.h

index 9c54a67..555d35c 100644 (file)
@@ -154,10 +154,6 @@ public:
        /** Clear some status_flags. */
        const basic & clearflag(unsigned f) const {flags &= ~f; return *this;}
 
-       /** Get relative precedence level (useful for implementing pretty-printed
-        *  output). */
-       unsigned get_precedence(void) const {return precedence;}
-
 protected:
        void ensure_if_modifiable(void) const;
        
index a903881..bd0c896 100644 (file)
@@ -214,7 +214,6 @@ ex clifford::simplify_ncmul(const exvector & v) const
 {
        exvector s;
        s.reserve(v.size());
-       unsigned rl = ex_to_clifford(v[0]).get_representation_label();
 
        // Remove superfluous ONEs
        exvector::const_iterator cit = v.begin(), citend = v.end();
@@ -266,7 +265,7 @@ ex clifford::simplify_ncmul(const exvector & v) const
                                const ex & ib = b.op(1);
                                if (ia.is_equal(ib)) {
                                        a = lorentz_g(ia, ib);
-                                       b = dirac_ONE(rl);
+                                       b = dirac_ONE(representation_label);
                                        something_changed = true;
                                }
                        }
@@ -275,7 +274,7 @@ ex clifford::simplify_ncmul(const exvector & v) const
        }
 
        if (s.size() == 0)
-               return clifford(diracone(), rl) * sign;
+               return clifford(diracone(), representation_label) * sign;
        if (something_changed)
                return nonsimplified_ncmul(s) * sign;
        else
@@ -478,25 +477,21 @@ ex canonicalize_clifford(const ex & e)
                        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++;
+               exvector::iterator it = v.begin(), next_to_last = v.end() - 1;
+               if (is_ex_of_type(it->op(0), diracgamma5))
+                       it++;
+               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);
                        }
-                       next_to_last--;
+                       it++;
                }
                return ncmul(v);
        }
index 9e719ec..f06805e 100644 (file)
@@ -156,7 +156,6 @@ ex color::simplify_ncmul(const exvector & v) const
 {
        exvector s;
        s.reserve(v.size());
-       unsigned rl = ex_to_color(v[0]).get_representation_label();
 
        // Remove superfluous ONEs
        exvector::const_iterator it = v.begin(), itend = v.end();
@@ -167,7 +166,7 @@ ex color::simplify_ncmul(const exvector & v) const
        }
 
        if (s.size() == 0)
-               return color(su3one(), rl);
+               return color(su3one(), representation_label);
        else
                return simplified_ncmul(s);
 }
@@ -357,8 +356,12 @@ bool su3d::contract_with(exvector::iterator self, exvector::iterator other, exve
        if (is_ex_exactly_of_type(other->op(0), su3d)) {
 
                // Find the dummy indices of the contraction
-               exvector dummy_indices;
-               dummy_indices = ex_to_indexed(*self).get_dummy_indices(ex_to_indexed(*other));
+               exvector self_indices = ex_to_indexed(*self).get_indices();
+               exvector other_indices = ex_to_indexed(*other).get_indices();
+               exvector all_indices = self_indices;
+               all_indices.insert(all_indices.end(), other_indices.begin(), other_indices.end());
+               exvector free_indices, dummy_indices;
+               find_free_and_dummy(all_indices, free_indices, dummy_indices);
 
                // d.abc*d.abc=40/3
                if (dummy_indices.size() == 3) {
@@ -368,11 +371,12 @@ bool su3d::contract_with(exvector::iterator self, exvector::iterator other, exve
 
                // 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]);
+                       exvector a;
+                       back_insert_iterator<exvector> ita(a);
+                       ita = set_difference(self_indices.begin(), self_indices.end(), dummy_indices.begin(), dummy_indices.end(), ita, ex_is_less());
+                       ita = set_difference(other_indices.begin(), other_indices.end(), dummy_indices.begin(), dummy_indices.end(), ita, ex_is_less());
+                       GINAC_ASSERT(a.size() == 2);
+                       *self = numeric(5, 3) * delta_tensor(a[0], a[1]);
                        *other = _ex1();
                        return true;
                }
index 90ac3a4..d54f519 100644 (file)
@@ -21,6 +21,7 @@
  */
 
 #include <stdexcept>
+#include <algorithm>
 
 #include "idx.h"
 #include "symbol.h"
@@ -443,29 +444,36 @@ bool is_dummy_pair(const ex & e1, const ex & 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)
+// Shaker sort is sufficient for the expected small number of indices
+template <class It, class Cmp>
+inline void shaker_sort(It first, It last, Cmp comp)
 {
-       // Nothing to sort if less than 2 elements
-       if (v.size() < 2)
+       if (first == last)
                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++;
+       --last;
+       if (first == last)
+               return;
+       It flag = first;
+       do {
+               It i;
+               for (i=last; i>first; --i) {
+                       if (comp(*i, i[-1])) {
+                               iter_swap(i-1, i);
+                               flag = i - 1;
+                       }
                }
-               it1++;
-       }
+               ++flag;
+               first = flag;
+               for (i=first; i<last; ++i) {
+                       if (comp(i[1], *i)) {
+                               iter_swap(i, i+1);
+                               flag = i + 1;
+                       }
+               }
+               last = flag - 1;
+       } while (first <= last);
 }
 
-
 void find_free_and_dummy(exvector::const_iterator it, exvector::const_iterator itend, exvector & out_free, exvector & out_dummy)
 {
        out_free.clear();
@@ -485,7 +493,7 @@ void find_free_and_dummy(exvector::const_iterator it, exvector::const_iterator i
        // 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);
+       shaker_sort(v.begin(), v.end(), ex_is_less());
 
        // Find dummy pairs and free indices
        it = v.begin(); itend = v.end();
@@ -506,27 +514,4 @@ void find_free_and_dummy(exvector::const_iterator it, exvector::const_iterator i
                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
index 07ff839..f8c1f06 100644 (file)
@@ -241,10 +241,6 @@ inline unsigned count_free_indices(const exvector & v)
        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
 
 #endif // ndef __GINAC_IDX_H__
index 435d952..f8d0903 100644 (file)
@@ -225,6 +225,12 @@ bool indexed::info(unsigned inf) const
        return inherited::info(inf);
 }
 
+struct idx_is_not : public binary_function<ex, unsigned, bool> {
+       bool operator() (const ex & e, unsigned inf) const {
+               return !(ex_to_idx(e).get_value().info(inf));
+       }
+};
+
 bool indexed::all_index_values_are(unsigned inf) const
 {
        // No indices? Then no property can be fulfilled
@@ -232,14 +238,7 @@ bool indexed::all_index_values_are(unsigned inf) const
                return false;
 
        // Check all indices
-       exvector::const_iterator it = seq.begin() + 1, itend = seq.end();
-       while (it != itend) {
-               GINAC_ASSERT(is_ex_of_type(*it, idx));
-               if (!ex_to_idx(*it).get_value().info(inf))
-                       return false;
-               it++;
-       }
-       return true;
+       return find_if(seq.begin() + 1, seq.end(), bind2nd(idx_is_not(), inf)) == seq.end();
 }
 
 int indexed::compare_same_type(const basic & other) const
@@ -454,15 +453,7 @@ static bool indices_consistent(const exvector & v1, const exvector & v2)
        if (v1.size() != v2.size())
                return false;
 
-       // And also the indices themselves
-       exvector::const_iterator ait = v1.begin(), aitend = v1.end(),
-                                bit = v2.begin(), bitend = v2.end();
-       while (ait != aitend) {
-               if (!ait->is_equal(*bit))
-                       return false;
-               ait++; bit++;
-       }
-       return true;
+       return equal(v1.begin(), v1.end(), v2.begin(), ex_is_equal());
 }
 
 exvector indexed::get_indices(void) const
@@ -546,14 +537,6 @@ 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
@@ -571,8 +554,6 @@ static ex rename_dummy_indices(const ex & e, exvector & global_dummy_indices, ex
        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
@@ -580,18 +561,13 @@ static ex rename_dummy_indices(const ex & e, exvector & global_dummy_indices, ex
                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++;
+                       if (find_if(global_dummy_indices.begin(), global_dummy_indices.end(), bind2nd(ex_is_equal(), *it)) == global_dummy_indices.end()) {
+                               global_dummy_indices.push_back(*it);
+                               global_size++;
+                               remaining--;
                        }
-                       global_dummy_indices.push_back(*it);
-                       global_size++;
-                       remaining--;
-found:         it++;
+                       it++;
                }
-               sort(global_dummy_indices.begin(), global_dummy_indices.end(), ex_is_less());
        }
 
        // Replace index symbols in expression
@@ -601,10 +577,11 @@ found:            it++;
        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))
+               if (!loc_sym.is_equal(glob_sym)) {
                        all_equal = false;
-               local_syms.append(loc_sym);
-               global_syms.append(glob_sym);
+                       local_syms.append(loc_sym);
+                       global_syms.append(glob_sym);
+               }
        }
        if (all_equal)
                return e;
@@ -732,14 +709,23 @@ contraction_done:
        }
 
        // Find free indices (concatenate them all and call find_free_and_dummy())
-       exvector un, local_dummy_indices;
+       // and all dummy indices that appear
+       exvector un, individual_dummy_indices;
        it1 = v.begin(); itend = v.end();
        while (it1 != itend) {
-               exvector free_indices_of_factor = it1->get_free_indices();
+               exvector free_indices_of_factor;
+               if (is_ex_of_type(*it1, indexed)) {
+                       exvector dummy_indices_of_factor;
+                       find_free_and_dummy(ex_to_indexed(*it1).seq.begin() + 1, ex_to_indexed(*it1).seq.end(), free_indices_of_factor, dummy_indices_of_factor);
+                       individual_dummy_indices.insert(individual_dummy_indices.end(), dummy_indices_of_factor.begin(), dummy_indices_of_factor.end());
+               } else
+                       free_indices_of_factor = it1->get_free_indices();
                un.insert(un.end(), free_indices_of_factor.begin(), free_indices_of_factor.end());
                it1++;
        }
+       exvector local_dummy_indices;
        find_free_and_dummy(un, free_indices, local_dummy_indices);
+       local_dummy_indices.insert(local_dummy_indices.end(), individual_dummy_indices.begin(), individual_dummy_indices.end());
 
        ex r;
        if (something_changed)
@@ -765,7 +751,7 @@ ex simplify_indexed(const ex & e, exvector & free_indices, exvector & dummy_indi
        ex e_expanded = e.expand();
 
        // Simplification of single indexed object: just find the free indices
-       // (and perform dummy index renaming if 
+       // and perform dummy index renaming
        if (is_ex_of_type(e_expanded, indexed)) {
                const indexed &i = ex_to_indexed(e_expanded);
                exvector local_dummy_indices;
index 8831ada..53929bc 100644 (file)
@@ -91,7 +91,6 @@ protected:
        unsigned row;             ///< number of rows
        unsigned col;             ///< number of columns
        exvector m;               ///< representation (cols indexed first)
-       static unsigned precedence;
 };
 
 
index 660a9fa..b13665e 100644 (file)
@@ -673,9 +673,7 @@ ex power::expand_add(const add & a, int n) const
                        cout << "k_cum[" << i << "]=" << k_cum[i] << endl;
                        cout << "upper_limit[" << i << "]=" << upper_limit[i] << endl;
                }
-               for (exvector::const_iterator cit=term.begin(); cit!=term.end(); ++cit) {
-                       cout << *cit << endl;
-               }
+               for_each(term.begin(), term.end(), ostream_iterator<ex>(cout, "\n"));
                cout << "end term" << endl;
                */
                
index 125d099..8cab203 100644 (file)
@@ -81,6 +81,7 @@ public: \
        classname(const classname & other); \
        const classname & operator=(const classname & other); \
        basic * duplicate() const; \
+       unsigned get_precedence(void) const {return precedence;} \
 protected: \
        void copy(const classname & other); \
        void destroy(bool call_parent); \
index 93408d4..7acb3a9 100644 (file)
@@ -59,13 +59,6 @@ unsigned log2(unsigned n)
 }
 #endif
 
-/** Append one exvector to another */
-void append_exvector_to_exvector(exvector & dest, const exvector & source)
-{
-       dest.reserve(dest.size() + source.size());
-       dest.insert(dest.end(), source.begin(), source.end());
-}
-
 //////////
 // `construct on first use' chest of numbers
 //////////
index 6edea86..512c97c 100644 (file)
@@ -28,6 +28,7 @@
 
 #include <string>
 #include <stdexcept>
+#include <functional>
 #if defined(HAVE_SSTREAM)
 #include <sstream>
 #elif defined(HAVE_STRSTREAM)
@@ -146,7 +147,14 @@ int permutation_sign(std::vector<T> s)
        return sigma;
 }
 
-void append_exvector_to_exvector(exvector & dest, const exvector & source);
+/* Function objects for STL sort() etc. */
+struct ex_is_less : public binary_function<ex, ex, bool> {
+       bool operator() (const ex &lh, const ex &rh) const { return lh.compare(rh) < 0; }
+};
+
+struct ex_is_equal : public binary_function<ex, ex, bool> {
+       bool operator() (const ex &lh, const ex &rh) const { return lh.is_equal(rh); }
+};
 
 // Collection of `construct on first use' wrappers for safely avoiding
 // internal object replication without running into the `static