[GiNaC-devel] Bug(?) in reposition_dummy_indices: test case and
patch
Sheplyakov Alexei
varg at theor.jinr.ru
Tue Aug 29 19:39:18 CEST 2006
Hello!
This simple program
#include <iostream>
#include <stdexcept>
#include <ginac/ginac.h>
using namespace std;
using namespace GiNaC;
int main(int argc, char** argv)
{
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);
ex test_cycl = indexed(symbol("A"), sy_cycl(), ev2);
test_cycl = test_cycl.simplify_indexed();
if (test_cycl.get_free_indices().size() != 0)
{
std::cerr << e << std::endl;
throw std::logic_error("Inconsistent indices in expression");
}
return 0;
}
fails (both with GiNaC 1.4 CVS and 1.3.x). I believe that the reason
is bug in indexed.cpp:reposition_dummy_indices() and propose attached
patch to fix it. Note that patch *seems* to be correct, but IMHO it
is somewhat ugly and probably inefficient. Could anyone suggest a better
solution?
Best regards,
Alexei.
--
All science is either physics or stamp collecting.
-------------- next part --------------
Index: ginac/indexed.cpp
===================================================================
RCS file: /home/cvs/GiNaC/ginac/indexed.cpp,v
retrieving revision 1.104
diff -u -r1.104 indexed.cpp
--- ginac/indexed.cpp 28 Jul 2006 22:34:28 -0000 1.104
+++ ginac/indexed.cpp 29 Aug 2006 17:24:02 -0000
@@ -628,6 +628,36 @@
}
}
+/** Exchange variance of dummy indicies, helper class used
+ * by reposition_dummy_indices */
+class exchange_varidx_variance_fct
+{
+ private:
+ varidx vi;
+ public:
+ explicit exchange_varidx_variance_fct(const varidx& vi_) : vi(vi_) { }
+ inline void operator()(ex& arg) const
+ {
+ if (arg.is_equal(vi) || arg.is_equal(vi.toggle_variance()))
+ arg = ex_to<varidx>(arg).toggle_variance();
+ }
+};
+
+/** Turn contraviariant varidx into covariant one, helper class used
+ * by reposition_dummy_indices */
+class lower_dummy_varidx_fct
+{
+ private:
+ varidx vi;
+ public:
+ explicit lower_dummy_varidx_fct(const varidx& vi_): vi(vi_) { }
+ inline void operator()(ex & arg) const
+ {
+ if (arg.is_equal(vi))
+ arg = vi.toggle_variance();
+ }
+};
+
/** Raise/lower dummy indices in a single indexed objects to canonicalize their
* variance.
*
@@ -638,26 +668,41 @@
bool reposition_dummy_indices(ex & e, exvector & variant_dummy_indices, exvector & moved_indices)
{
bool something_changed = false;
-
+ 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) {
+ exvector::iterator it2, it2start, it2end;
+ for (it2start = seq.begin(), it2end = seq.end(), it2 = it2start + 1; it2 != it2end; ++it2)
+ {
if (!is_exactly_a<varidx>(*it2))
continue;
exvector::iterator vit, vitend;
- 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);
+ 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())
+ {
+ /*
+ * N.B. we can't 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
+ */
+ std::for_each(seq.begin(), seq.end(),
+ exchange_varidx_variance_fct(ex_to<varidx>(*it2)));
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();
+ it2 = seq.begin() + (it2 - it2start);
+ it2start = seq.begin();
+ it2end = seq.end();
}
moved_indices.push_back(*vit);
variant_dummy_indices.erase(vit);
@@ -668,11 +713,12 @@
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);
+ std::for_each(seq.begin(), seq.end(),
+ lower_dummy_varidx_fct(ex_to<varidx>(*it2)));
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();
+ it2 = seq.begin() + (it2 - it2start);
+ it2start = seq.begin();
+ it2end = seq.end();
}
goto next_index;
}
@@ -681,6 +727,13 @@
next_index: ;
}
+ if (something_changed)
+ {
+ indexed ei_ = ex_to<indexed>(e);
+ ei_.seq = seq;
+ e = ei_;
+ }
+
return something_changed;
}
Index: check/exam_paranoia.cpp
===================================================================
RCS file: /home/cvs/GiNaC/check/exam_paranoia.cpp,v
retrieving revision 1.20
diff -u -r1.20 exam_paranoia.cpp
--- check/exam_paranoia.cpp 1 May 2005 18:23:26 -0000 1.20
+++ check/exam_paranoia.cpp 29 Aug 2006 17:24:02 -0000
@@ -437,6 +437,33 @@
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 @@
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;
-------------- next part --------------
A non-text attachment was scrubbed...
Name: not available
Type: application/pgp-signature
Size: 827 bytes
Desc: Digital signature
Url : http://www.cebix.net/pipermail/ginac-devel/attachments/20060829/95547bca/attachment.pgp
More information about the GiNaC-devel
mailing list