Hi, Chris! On Fri, Oct 20, 2006 at 10:51:21AM +0200, Chris Dams wrote:
Yes. However, if we have a traceless tensor both T~mu.mu and T.mu~mu should evaluated themselves into 0. It does not make sense to only have one of them evaluate into zero. The point I was trying to make is that tensors with some index lowered/raised always (?) have the same properties so that an evaluate step that aplies to the version with an up-index should also aply to the one with a down-index.
I think it *does* make sense to replace only one of combinations and rely on simplify_indexed to get indices canonicalized. Otherwise one need to duplicate quite a large portion of code from indexed.cpp and symmetry.cpp (to decide what combinations are equivalent). [snipped]
Oh well, the thing seemed rather complicated to me and didn't do all simplifications possible, as noted in http://www.ginac.de/pipermail/ginac-devel/2006-August/001055.html. I now more or less applied it, but with some modifications. I hope that I haven't broken anything ;-). (1) Now that the raising and lowering of indices is done on an exvector, we do not need to raise/lower both indices at the same time. Another reason why this isn't necessary anymore is that indices that occur in dummy pairs within the indexed object e should already have been removed from the vector variant_dummy_indices. I think this means that the for_each that you had can be omitted.
I agree.
(2) It is not correct to return the indexed object by doing
indexed ei_ = ex_to<indexed>(e); ei_.seq = seq; e = ei_;
This should be e = ex_to<indexed>(e).thiscontainer(seq); because we may be dealing with something that inherits from indexed rather than an indexed object.
I didn't like that chunk myself exactly for this reason, but I didn't know the correct way.
I also simplified the clifford exam by removing a huge macro.
A patch is in CVS. Works fine for me. I've re-diffed it for 1.3 branch, could you please apply it?
Best regards, Alexei. [PATCH] reposition_dummy_indices: fix bugs w.r.t. raising/lowering dummy indices. Now reposition_dummy_indices works fine if the argument has non-trivial symmetry properties and/or evaluation method (which triggers re-evaluation in course of exchanging of dummy indices). Thus, simplify_indexed gives correct results for such expressions. Backported from the main branch. --- check/exam_paranoia.cpp | 28 ++++++++++++++++ ginac/indexed.cpp | 82 ++++++++++++++++++++++++++++++++++++++++------- 2 files changed, 97 insertions(+), 13 deletions(-) diff --git a/check/exam_paranoia.cpp b/check/exam_paranoia.cpp index 7c82e48..0ab6491 100644 --- a/check/exam_paranoia.cpp +++ b/check/exam_paranoia.cpp @@ -437,6 +437,33 @@ static unsigned exam_paranoia16() return result; } +// Bug in reposition_dummy_indices() could result in correct expression +// turned into one with inconsistent indices. Fixed on Aug 29, 2006 +static unsigned exam_paranoia17() +{ + varidx mu1(symbol("mu1"), 4); + varidx mu2(symbol("mu2"), 4); + varidx mu3(symbol("mu3"), 4); + varidx mu4(symbol("mu4"), 4); + varidx mu5(symbol("mu5"), 4); + varidx mu6(symbol("mu6"), 4); + + exvector ev2; + ev2.push_back(mu3.toggle_variance()); + ev2.push_back(mu6); + ev2.push_back(mu5.toggle_variance()); + ev2.push_back(mu6.toggle_variance()); + ev2.push_back(mu5); + ev2.push_back(mu3); + // notice: all indices are contracted ... + + ex test_cycl = indexed(symbol("A"), sy_cycl(), ev2); + test_cycl = test_cycl.simplify_indexed(); + // ... so there should be zero free indices in the end. + return test_cycl.get_free_indices().size(); +} + + unsigned exam_paranoia() { unsigned result = 0; @@ -460,6 +487,7 @@ unsigned exam_paranoia() result += exam_paranoia14(); cout << '.' << flush; result += exam_paranoia15(); cout << '.' << flush; result += exam_paranoia16(); cout << '.' << flush; + result += exam_paranoia17(); cout << '.' << flush; if (!result) { cout << " passed " << endl; diff --git a/ginac/indexed.cpp b/ginac/indexed.cpp index 05c3b90..2de726b 100644 --- a/ginac/indexed.cpp +++ b/ginac/indexed.cpp @@ -632,10 +632,60 @@ bool reposition_dummy_indices(ex & e, ex { bool something_changed = false; + // Find dummy symbols that occur twice in the same indexed object. + exvector local_var_dummies; + local_var_dummies.reserve(e.nops()/2); + for (size_t i=1; i<e.nops(); ++i) { + if (!is_a<varidx>(e.op(i))) + continue; + for (size_t j=i+1; j<e.nops(); ++j) { + if (is_dummy_pair(e.op(i), e.op(j))) { + local_var_dummies.push_back(e.op(i)); + for (exvector::iterator k = variant_dummy_indices.begin(); + k!=variant_dummy_indices.end(); ++k) { + if (e.op(i).op(0) == k->op(0)) { + variant_dummy_indices.erase(k); + break; + } + } + break; + } + } + } + + // In the case where a dummy symbol occurs twice in the same indexed object + // we try all posibilities of raising/lowering and keep the least one in + // the sense of ex_is_less. + ex optimal_e = e; + size_t numpossibs = 1 << local_var_dummies.size(); + for (size_t i=0; i<numpossibs; ++i) { + ex try_e = e; + for (size_t j=0; j<local_var_dummies.size(); ++j) { + exmap m; + if (1<<j & i) { + ex curr_idx = local_var_dummies[j]; + ex curr_toggle = ex_to<varidx>(curr_idx).toggle_variance(); + m[curr_idx] = curr_toggle; + m[curr_toggle] = curr_idx; + } + try_e = e.subs(m, subs_options::no_pattern); + } + if(ex_is_less()(try_e, optimal_e)) + { optimal_e = try_e; + something_changed = true; + } + } + e = optimal_e; + + if (!is_a<indexed>(e)) + return true; + + exvector seq = ex_to<indexed>(e).seq; + // If a dummy index is encountered for the first time in the // product, pull it up, otherwise, pull it down - exvector::const_iterator it2, it2start, it2end; - for (it2start = ex_to<indexed>(e).seq.begin(), it2end = ex_to<indexed>(e).seq.end(), it2 = it2start + 1; it2 != it2end; ++it2) { + for (exvector::iterator it2 = seq.begin()+1, it2end = seq.end(); + it2 != it2end; ++it2) { if (!is_exactly_a<varidx>(*it2)) continue; @@ -643,14 +693,20 @@ bool reposition_dummy_indices(ex & e, ex for (vit = variant_dummy_indices.begin(), vitend = variant_dummy_indices.end(); vit != vitend; ++vit) { if (it2->op(0).is_equal(vit->op(0))) { if (ex_to<varidx>(*it2).is_covariant()) { - e = e.subs(lst( - *it2 == ex_to<varidx>(*it2).toggle_variance(), - ex_to<varidx>(*it2).toggle_variance() == *it2 - ), subs_options::no_pattern); + /* + * N.B. we don't want to use + * + * e = e.subs(lst( + * *it2 == ex_to<varidx>(*it2).toggle_variance(), + * ex_to<varidx>(*it2).toggle_variance() == *it2 + * ), subs_options::no_pattern); + * + * since this can trigger non-trivial repositioning of indices, + * e.g. due to non-trivial symmetry properties of e, thus + * invalidating iterators + */ + *it2 = ex_to<varidx>(*it2).toggle_variance(); something_changed = true; - it2 = ex_to<indexed>(e).seq.begin() + (it2 - it2start); - it2start = ex_to<indexed>(e).seq.begin(); - it2end = ex_to<indexed>(e).seq.end(); } moved_indices.push_back(*vit); variant_dummy_indices.erase(vit); @@ -661,11 +717,8 @@ bool reposition_dummy_indices(ex & e, ex for (vit = moved_indices.begin(), vitend = moved_indices.end(); vit != vitend; ++vit) { if (it2->op(0).is_equal(vit->op(0))) { if (ex_to<varidx>(*it2).is_contravariant()) { - e = e.subs(*it2 == ex_to<varidx>(*it2).toggle_variance(), subs_options::no_pattern); + *it2 = ex_to<varidx>(*it2).toggle_variance(); something_changed = true; - it2 = ex_to<indexed>(e).seq.begin() + (it2 - it2start); - it2start = ex_to<indexed>(e).seq.begin(); - it2end = ex_to<indexed>(e).seq.end(); } goto next_index; } @@ -674,6 +727,9 @@ bool reposition_dummy_indices(ex & e, ex next_index: ; } + if (something_changed) + e = ex_to<indexed>(e).thiscontainer(seq); + return something_changed; } -- 1.4.2.3 -- All science is either physics or stamp collecting.