From 2967504b9188e7169bf0d947089ddf879140614d Mon Sep 17 00:00:00 2001 From: Christian Bauer Date: Tue, 21 Jan 2003 22:36:58 +0000 Subject: [PATCH] fixed bugs in dummy index symmetrization and clifford contraction --- ginac/clifford.cpp | 5 +- ginac/indexed.cpp | 156 +++++++++++++++++++++++++++++++++------------ 2 files changed, 120 insertions(+), 41 deletions(-) diff --git a/ginac/clifford.cpp b/ginac/clifford.cpp index 5a386e7d..262becb6 100644 --- a/ginac/clifford.cpp +++ b/ginac/clifford.cpp @@ -209,7 +209,10 @@ bool diracgamma::contract_with(exvector::iterator self, exvector::iterator other GINAC_ASSERT(is_a(*other)); GINAC_ASSERT(is_a(self->op(0))); unsigned char rl = ex_to(*self).get_representation_label(); + ex dim = ex_to(self->op(1)).get_dim(); + if (other->nops() > 1) + dim = minimal_dim(dim, ex_to(self->op(1)).get_dim()); if (is_a(*other)) { @@ -440,7 +443,7 @@ ex clifford::eval_ncmul(const exvector & v) const } else if (!a_is_diracgamma && !b_is_diracgamma && ag.is_equal(bg)) { // a\ a\ -> a^2 - varidx ix((new symbol)->setflag(status_flags::dynallocated), ex_to(a.op(1)).get_dim()); + varidx ix((new symbol)->setflag(status_flags::dynallocated), ex_to(a.op(1)).minimal_dim(ex_to(b.op(1)))); a = indexed(ag, ix) * indexed(ag, ix.toggle_variance()); b = dirac_ONE(representation_label); something_changed = true; diff --git a/ginac/indexed.cpp b/ginac/indexed.cpp index ae2328ad..366db44d 100644 --- a/ginac/indexed.cpp +++ b/ginac/indexed.cpp @@ -21,6 +21,7 @@ */ #include +#include #include #include "indexed.h" @@ -388,6 +389,22 @@ ex indexed::derivative(const symbol & s) const // global functions ////////// +struct idx_is_equal_ignore_dim : public std::binary_function { + bool operator() (const ex &lh, const ex &rh) const + { + if (lh.is_equal(rh)) + return true; + else + try { + // Replacing the dimension might cause an error (e.g. with + // index classes that only work in a fixed number of dimensions) + return lh.is_equal(ex_to(rh).replace_dim(ex_to(lh).get_dim())); + } catch (...) { + return false; + } + } +}; + /** Check whether two sorted index vectors are consistent (i.e. equal). */ static bool indices_consistent(const exvector & v1, const exvector & v2) { @@ -395,7 +412,7 @@ static bool indices_consistent(const exvector & v1, const exvector & v2) if (v1.size() != v2.size()) return false; - return equal(v1.begin(), v1.end(), v2.begin(), ex_is_equal()); + return equal(v1.begin(), v1.end(), v2.begin(), idx_is_equal_ignore_dim()); } exvector indexed::get_indices(void) const @@ -833,17 +850,11 @@ contraction_done: * obtained during the simplification of sums. */ class symminfo { public: - symminfo() {} + symminfo() : num(0) {} ~symminfo() {} - symminfo(const ex & symmterm_, const ex & orig_) + symminfo(const ex & symmterm_, const ex & orig_, unsigned num_) : orig(orig_), num(num_) { - if (is_exactly_a(orig_) && is_exactly_a(orig_.op(orig_.nops()-1))) { - ex tmp = orig_.op(orig_.nops()-1); - orig = orig_ / tmp; - } else - orig = orig_; - if (is_exactly_a(symmterm_) && is_exactly_a(symmterm_.op(symmterm_.nops()-1))) { coeff = symmterm_.op(symmterm_.nops()-1); symmterm = symmterm_ / coeff; @@ -853,40 +864,44 @@ public: } } - symminfo(const symminfo & other) - { - symmterm = other.symmterm; - coeff = other.coeff; - orig = other.orig; - } + ex symmterm; /**< symmetrized term */ + ex coeff; /**< coefficient of symmetrized term */ + ex orig; /**< original term */ + unsigned num; /**< how many symmetrized terms resulted from the original term */ +}; - const symminfo & operator=(const symminfo & other) +class symminfo_is_less_by_symmterm { +public: + bool operator() (const symminfo & si1, const symminfo & si2) const { - if (this != &other) { - symmterm = other.symmterm; - coeff = other.coeff; - orig = other.orig; - } - return *this; + int comp = si1.symmterm.compare(si2.symmterm); + if (comp < 0) return true; +#if 0 + if (comp > 0) return false; + comp = si1.orig.compare(si2.orig); + if (comp < 0) return true; + if (comp > 0) return false; + comp = si1.coeff.compare(si2.coeff); + if (comp < 0) return true; +#endif + return false; } - - ex symmterm; - ex coeff; - ex orig; }; -class symminfo_is_less { +class symminfo_is_less_by_orig { public: - bool operator() (const symminfo & si1, const symminfo & si2) + bool operator() (const symminfo & si1, const symminfo & si2) const { - int comp = si1.symmterm.compare(si2.symmterm); + int comp = si1.orig.compare(si2.orig); if (comp < 0) return true; +#if 0 if (comp > 0) return false; - comp = si1.orig.compare(si2.orig); + comp = si1.symmterm.compare(si2.symmterm); if (comp < 0) return true; if (comp > 0) return false; comp = si1.coeff.compare(si2.coeff); if (comp < 0) return true; +#endif return false; } }; @@ -938,8 +953,12 @@ ex simplify_indexed(const ex & e, exvector & free_indices, exvector & dummy_indi sum = term; first = false; } else { - if (!indices_consistent(free_indices, free_indices_of_term)) - throw (std::runtime_error("simplify_indexed: inconsistent indices in sum")); + if (!indices_consistent(free_indices, free_indices_of_term)) { + std::ostringstream s; + s << "simplify_indexed: inconsistent indices in sum: "; + s << exprseq(free_indices) << " vs. " << exprseq(free_indices_of_term); + throw (std::runtime_error(s.str())); + } if (is_a(sum) && is_a(term)) sum = ex_to(sum.op(0)).add_indexed(sum, term); else @@ -968,24 +987,81 @@ ex simplify_indexed(const ex & e, exvector & free_indices, exvector & dummy_indi std::vector v; for (int i=0; i(sum_symm)) for (int j=0; j::iterator i=v.begin(); i!=v.end(); ) { - std::vector::iterator j = i; - for (j++; j!=v.end() && i->symmterm == j->symmterm; j++) ; - for (std::vector::iterator k=i; k!=j; k++) - result.push_back((k->coeff)*(i->orig)); + std::sort(v.begin(), v.end(), symminfo_is_less_by_symmterm()); + + std::vector v_pass2; + for (std::vector::const_iterator i=v.begin(); i!=v.end(); ) { + + // Combine equal terms + std::vector::const_iterator j = i + 1; + if (j != v.end() && j->symmterm == i->symmterm) { + + // More than one term, collect the coefficients + ex coeff = i->coeff; + while (j != v.end() && j->symmterm == i->symmterm) { + coeff += j->coeff; + j++; + } + + // Add combined term to result + if (!coeff.is_zero()) + result.push_back(coeff * i->symmterm); + + } else { + + // Single term, store for second pass + v_pass2.push_back(*i); + } + i = j; } + + // Were there any remaining terms that didn't get combined? + if (v_pass2.size() > 0) { + + // Yes, sort them by their original term + std::sort(v_pass2.begin(), v_pass2.end(), symminfo_is_less_by_orig()); + + for (std::vector::const_iterator i=v_pass2.begin(); i!=v_pass2.end(); ) { + + // How many symmetrized terms of this original term are left? + unsigned num = 1; + std::vector::const_iterator j = i + 1; + while (j != v_pass2.end() && j->orig == i->orig) { + num++; + j++; + } + + if (num == i->num) { + + // All terms left, then add the original term to the result + result.push_back(i->orig); + + } else { + + // Some terms were combined with others, add up the remaining symmetrized terms + std::vector::const_iterator k; + for (k=i; k!=j; k++) + result.push_back(k->coeff * k->symmterm); + } + + i = j; + } + } + + // Add all resulting terms ex sum_symm = (new add(result))->setflag(status_flags::dynallocated); if (sum_symm.is_zero()) free_indices.clear(); -- 2.49.0