From df7b9291027e0e5bda65e07fe251469ef964e704 Mon Sep 17 00:00:00 2001 From: Christian Bauer Date: Sat, 19 May 2001 00:42:07 +0000 Subject: [PATCH] - dummy index renamer didn't account for internal dummy indices of objects 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 | 4 --- ginac/clifford.cpp | 37 +++++++++++------------- ginac/color.cpp | 22 ++++++++------ ginac/idx.cpp | 69 +++++++++++++++++--------------------------- ginac/idx.h | 4 --- ginac/indexed.cpp | 72 +++++++++++++++++++--------------------------- ginac/matrix.h | 1 - ginac/power.cpp | 4 +-- ginac/registrar.h | 1 + ginac/utils.cpp | 7 ----- ginac/utils.h | 10 ++++++- 11 files changed, 96 insertions(+), 135 deletions(-) diff --git a/ginac/basic.h b/ginac/basic.h index 9c54a67e..555d35cf 100644 --- a/ginac/basic.h +++ b/ginac/basic.h @@ -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; diff --git a/ginac/clifford.cpp b/ginac/clifford.cpp index a9038818..bd0c8961 100644 --- a/ginac/clifford.cpp +++ b/ginac/clifford.cpp @@ -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); } diff --git a/ginac/color.cpp b/ginac/color.cpp index 9e719ec0..f06805e9 100644 --- a/ginac/color.cpp +++ b/ginac/color.cpp @@ -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 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; } diff --git a/ginac/idx.cpp b/ginac/idx.cpp index 90ac3a42..d54f5192 100644 --- a/ginac/idx.cpp +++ b/ginac/idx.cpp @@ -21,6 +21,7 @@ */ #include +#include #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 +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; iis_equal(*bit)) { - found = true; - break; - } - bit++; - } - if (!found) - ret.push_back(*ait); - ait++; - } - - return ret; -} - } // namespace GiNaC diff --git a/ginac/idx.h b/ginac/idx.h index 07ff8392..f8c1f060 100644 --- a/ginac/idx.h +++ b/ginac/idx.h @@ -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__ diff --git a/ginac/indexed.cpp b/ginac/indexed.cpp index 435d9527..f8d09030 100644 --- a/ginac/indexed.cpp +++ b/ginac/indexed.cpp @@ -225,6 +225,12 @@ bool indexed::info(unsigned inf) const return inherited::info(inf); } +struct idx_is_not : public binary_function { + 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; iget_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; diff --git a/ginac/matrix.h b/ginac/matrix.h index 8831ada8..53929bc5 100644 --- a/ginac/matrix.h +++ b/ginac/matrix.h @@ -91,7 +91,6 @@ protected: unsigned row; ///< number of rows unsigned col; ///< number of columns exvector m; ///< representation (cols indexed first) - static unsigned precedence; }; diff --git a/ginac/power.cpp b/ginac/power.cpp index 660a9fab..b13665e1 100644 --- a/ginac/power.cpp +++ b/ginac/power.cpp @@ -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(cout, "\n")); cout << "end term" << endl; */ diff --git a/ginac/registrar.h b/ginac/registrar.h index 125d0998..8cab2031 100644 --- a/ginac/registrar.h +++ b/ginac/registrar.h @@ -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); \ diff --git a/ginac/utils.cpp b/ginac/utils.cpp index 93408d4c..7acb3a98 100644 --- a/ginac/utils.cpp +++ b/ginac/utils.cpp @@ -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 ////////// diff --git a/ginac/utils.h b/ginac/utils.h index 6edea866..512c97c3 100644 --- a/ginac/utils.h +++ b/ginac/utils.h @@ -28,6 +28,7 @@ #include #include +#include #if defined(HAVE_SSTREAM) #include #elif defined(HAVE_STRSTREAM) @@ -146,7 +147,14 @@ int permutation_sign(std::vector 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 { + bool operator() (const ex &lh, const ex &rh) const { return lh.compare(rh) < 0; } +}; + +struct ex_is_equal : public binary_function { + 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 -- 2.49.0