Bugs w.r.t. raising/lowering dummy indices when doing simplify_indexed.

diff --git a/check/exam_clifford.cpp b/check/exam_clifford.cpp
index d6544eda5a3461480600d9b31a908167df09c830..e8ee3f70939e2131ede1bc549c554f18eaabd5a2 100644 (file)
--- a/check/exam_clifford.cpp
+++ b/check/exam_clifford.cpp
@@ -64,9 +64,14 @@ static unsigned check_equal_simplify_term(const ex & e1, const ex & e2, idx & mu
        ex e = expand_dummy_sum(normal(simplify_indexed(e1) - e2), true);
 
        for (int j=0; j<4; j++) {
-               ex esub = e.subs(lst(is_a<varidx>(mu) ? 
-                                                        mu == idx(j, mu.get_dim()), ex_to<varidx>(mu).toggle_variance() == idx(j, mu.get_dim()) 
-                                                        : mu == idx(j, mu.get_dim())));
+               ex esub = e.subs(
+                               is_a<varidx>(mu)
+                                       ? lst (
+                                                       mu == idx(j, mu.get_dim()),
+                                                       ex_to<varidx>(mu).toggle_variance() == idx(j, mu.get_dim())
+                                               )
+                                       : lst(mu == idx(j, mu.get_dim()))
+                       );
                if (!(canonicalize_clifford(esub).is_zero())) {
                        clog << "simplify_indexed(" << e1 << ") - (" << e2 << ") erroneously returned "
                                 << canonicalize_clifford(esub) << " instead of 0 for mu=" << j << endl;
@@ -312,151 +317,139 @@ static unsigned clifford_check5()
  * the cases when used indexes have or have not variance.
  * To this end we recycle the code through the following macros */
 
-#define CHECK6(IDX,TOGGLE) {IDX v(symbol("v"), 4), nu(symbol("nu"), 4), mu(symbol("mu"), 4), \
-              psi(symbol("psi"),4), lam(symbol("lambda"), 4),\
-              xi(symbol("xi"), 4),  rho(symbol("rho"),4);\
-\
-/* checks general identities and contractions for clifford_unit*/\
-       e = dirac_ONE(2) * clifford_unit(mu, A, 2) * dirac_ONE(2);\
-       result += check_equal(e, clifford_unit(mu, A, 2));\
-\
-       e = clifford_unit(idx(2, 4), A) * clifford_unit(idx(1, 4), A)\
-         * clifford_unit(idx(1, 4), A) * clifford_unit(idx(2, 4), A);\
-       result += check_equal(e, A(1, 1) * A(2, 2) * dirac_ONE());\
-\
-       e = clifford_unit(IDX(2, 4), A) * clifford_unit(IDX(1, 4), A)\
-         * clifford_unit(IDX(1, 4), A) * clifford_unit(IDX(2, 4), A);\
-       result += check_equal(e, A(1, 1) * A(2, 2) * dirac_ONE());\
-\
-       e = clifford_unit(nu, A) * clifford_unit(nu TOGGLE, A);\
-       result += check_equal_simplify(e, A.trace() * dirac_ONE());\
-\
-       e = clifford_unit(nu, A) * clifford_unit(nu, A);\
-       result += check_equal_simplify(e, indexed(A_symm, sy_symm(), nu, nu) * dirac_ONE());\
-\
-       e = clifford_unit(nu, A) * clifford_unit(nu TOGGLE, A) * clifford_unit(mu, A);\
-       result += check_equal_simplify(e, A.trace() * clifford_unit(mu, A));\
-\
-       e = clifford_unit(nu, A) * clifford_unit(mu, A) * clifford_unit(nu TOGGLE, A);\
-       \
-       result += check_equal_simplify_term(e,  2 * indexed(A_symm, sy_symm(), nu TOGGLE, mu) *clifford_unit(nu, A)-A.trace()*clifford_unit(mu, A), mu);\
-\
-       e = clifford_unit(nu, A) * clifford_unit(nu TOGGLE, A)\
-         * clifford_unit(mu, A) * clifford_unit(mu TOGGLE, A);\
-       result += check_equal_simplify(e, pow(A.trace(), 2) * dirac_ONE());\
-\
-       e = clifford_unit(mu, A) * clifford_unit(nu, A)\
-         * clifford_unit(nu TOGGLE, A) * clifford_unit(mu TOGGLE, A);\
-       result += check_equal_simplify(e, pow(A.trace(), 2)  * dirac_ONE());\
-\
-       e = clifford_unit(mu, A) * clifford_unit(nu, A)\
-         * clifford_unit(mu TOGGLE, A) * clifford_unit(nu TOGGLE, A);\
-\
-       result += check_equal_simplify_term2(e, 2*indexed(A_symm, sy_symm(), nu TOGGLE, mu TOGGLE) * clifford_unit(mu, A) * clifford_unit(nu, A) - pow(A.trace(), 2)*dirac_ONE());\
-\
-       e = clifford_unit(mu TOGGLE, A) * clifford_unit(nu, A)\
-         * clifford_unit(mu, A) * clifford_unit(nu TOGGLE, A);\
-\
-       result += check_equal_simplify_term2(e, 2*indexed(A_symm, nu, mu) * clifford_unit(mu TOGGLE, A) * clifford_unit(nu TOGGLE, A) - pow(A.trace(), 2)*dirac_ONE());\
-\
-       e = clifford_unit(nu TOGGLE, A) * clifford_unit(rho TOGGLE, A)\
-         * clifford_unit(mu, A) * clifford_unit(rho, A) * clifford_unit(nu, A);\
-       e = e.simplify_indexed().collect(clifford_unit(mu, A));\
-       \
-       result += check_equal_simplify_term(e, 4* indexed(A_symm, sy_symm(), nu TOGGLE,  rho)*indexed(A_symm, sy_symm(), rho TOGGLE, mu) *clifford_unit(nu, A) \
-                                           - 2*A.trace() * (clifford_unit(rho, A) * indexed(A_symm, sy_symm(), rho TOGGLE, mu) \
-                                                            + clifford_unit(nu, A) * indexed(A_symm, sy_symm(), nu TOGGLE, mu)) + pow(A.trace(),2)* clifford_unit(mu, A), mu);\
-\
-       e = clifford_unit(nu TOGGLE, A) * clifford_unit(rho, A)\
-         * clifford_unit(mu, A) * clifford_unit(rho TOGGLE, A) * clifford_unit(nu, A);\
-       e = e.simplify_indexed().collect(clifford_unit(mu, A));\
-       \
-       result += check_equal_simplify_term(e, 4* indexed(A_symm, sy_symm(), nu TOGGLE,  rho)*indexed(A_symm, sy_symm(), rho TOGGLE, mu) *clifford_unit(nu, A) \
-                                           - 2*A.trace() * (clifford_unit(rho, A) * indexed(A_symm, sy_symm(), rho TOGGLE, mu) \
-                                                            + clifford_unit(nu, A) * indexed(A_symm, sy_symm(), nu TOGGLE, mu)) + pow(A.trace(),2)* clifford_unit(mu, A), mu);\
-\
-       e = clifford_unit(mu, A) * clifford_unit(nu, A) + clifford_unit(nu, A) * clifford_unit(mu, A);\
-       result += check_equal(canonicalize_clifford(e), 2*dirac_ONE()*indexed(A_symm, sy_symm(), mu, nu));\
-\
-       e = (clifford_unit(mu, A) * clifford_unit(nu, A) * clifford_unit(lam, A)\
-                + clifford_unit(nu, A) * clifford_unit(lam, A) * clifford_unit(mu, A)\
-                + clifford_unit(lam, A) * clifford_unit(mu, A) * clifford_unit(nu, A)\
-                - clifford_unit(nu, A) * clifford_unit(mu, A) * clifford_unit(lam, A)\
-                - clifford_unit(lam, A) * clifford_unit(nu, A) * clifford_unit(mu, A)\
-                - clifford_unit(mu, A) * clifford_unit(lam, A) * clifford_unit(nu, A)) / 6\
-               + indexed(A_symm, sy_symm(), mu, nu) * clifford_unit(lam, A)\
-               - indexed(A_symm, sy_symm(), mu, lam) * clifford_unit(nu, A)\
-               + indexed(A_symm, sy_symm(), nu, lam) * clifford_unit(mu, A)\
-               - clifford_unit(mu, A) * clifford_unit(nu, A) * clifford_unit(lam, A);\
-       result += check_equal(canonicalize_clifford(e), 0);\
-\
-/* lst_to_clifford() and clifford_inverse()  check*/\
-       realsymbol x("x"), y("y"), t("t"), z("z");\
-       \
-       ex c = clifford_unit(nu, A, 1);\
-       e = lst_to_clifford(lst(t, x, y, z), mu, A, 1) * lst_to_clifford(lst(1, 2, 3, 4), c);\
-       e1 = clifford_inverse(e);\
-       result += check_equal_simplify_term2((e*e1).simplify_indexed(), dirac_ONE(1));\
-\
-/* Moebius map (both forms) checks for symmetric metrics only */\
-       matrix M1(2, 2),  M2(2, 2);\
-       c = clifford_unit(nu, A);\
-       \
-       e = clifford_moebius_map(0, dirac_ONE(), \
-                                                        dirac_ONE(), 0, lst(t, x, y, z), A); \
-/* this is just the inversion*/\
-       M1 = 0, dirac_ONE(),\
-               dirac_ONE(), 0;\
-       e1 = clifford_moebius_map(M1, lst(t, x, y, z), A); \
-/* the inversion again*/\
-       result += check_equal_lst(e, e1);\
-       \
-       e1 = clifford_to_lst(clifford_inverse(lst_to_clifford(lst(t, x, y, z), mu, A)), c);\
-       result += check_equal_lst(e, e1);\
-       \
-       e = clifford_moebius_map(dirac_ONE(), lst_to_clifford(lst(1, 2, 3, 4), nu, A), \
-                                                        0, dirac_ONE(), lst(t, x, y, z), A); \
-/*this is just a shift*/\
-       M2 = dirac_ONE(), lst_to_clifford(lst(1, 2, 3, 4), c),\
-               0, dirac_ONE();\
-       e1 = clifford_moebius_map(M2, lst(t, x, y, z), c); \
-/* the same shift*/\
-       result += check_equal_lst(e, e1);\
-               \
-       result += check_equal(e, lst(t+1, x+2, y+3, z+4));\
-       \
-/* Check the group law for Moebius maps */\
-       e = clifford_moebius_map(M1, ex_to<lst>(e1), c); \
-/*composition of M1 and M2*/\
-       e1 = clifford_moebius_map(M1.mul(M2), lst(t, x, y, z), c); \
-/* the product M1*M2*/\
-       result += check_equal_lst(e, e1);}
-
-static unsigned clifford_check6(const matrix & A)
+template <typename IDX> unsigned clifford_check6(const matrix &A)
 {
+       unsigned result = 0;
+
        matrix A_symm(4,4), A2(4, 4);
        A_symm = A.add(A.transpose()).mul(half);
        A2 = A_symm.mul(A_symm);
-       
+
+       IDX v(symbol("v"), 4), nu(symbol("nu"), 4), mu(symbol("mu"), 4),
+              psi(symbol("psi"),4), lam(symbol("lambda"), 4),
+              xi(symbol("xi"), 4),  rho(symbol("rho"),4);
+       ex mu_TOGGLE = is_a<varidx>(mu) ? ex_to<varidx>(mu).toggle_variance() : mu;
+       ex nu_TOGGLE = is_a<varidx>(nu) ? ex_to<varidx>(nu).toggle_variance() : nu;
+       ex rho_TOGGLE
+               = is_a<varidx>(rho) ? ex_to<varidx>(rho).toggle_variance() : rho;
+
        ex e, e1;
-       int result = 0;
+
+/* checks general identities and contractions for clifford_unit*/
+       e = dirac_ONE(2) * clifford_unit(mu, A, 2) * dirac_ONE(2);
+       result += check_equal(e, clifford_unit(mu, A, 2));
+
+       e = clifford_unit(idx(2, 4), A) * clifford_unit(idx(1, 4), A)
+         * clifford_unit(idx(1, 4), A) * clifford_unit(idx(2, 4), A);
+       result += check_equal(e, A(1, 1) * A(2, 2) * dirac_ONE());
+
+       e = clifford_unit(IDX(2, 4), A) * clifford_unit(IDX(1, 4), A)
+         * clifford_unit(IDX(1, 4), A) * clifford_unit(IDX(2, 4), A);
+       result += check_equal(e, A(1, 1) * A(2, 2) * dirac_ONE());
+
+       e = clifford_unit(nu, A) * clifford_unit(nu_TOGGLE, A);
+       result += check_equal_simplify(e, A.trace() * dirac_ONE());
+
+       e = clifford_unit(nu, A) * clifford_unit(nu, A);
+       result += check_equal_simplify(e, indexed(A_symm, sy_symm(), nu, nu) * dirac_ONE());
+
+       e = clifford_unit(nu, A) * clifford_unit(nu_TOGGLE, A) * clifford_unit(mu, A);
+       result += check_equal_simplify(e, A.trace() * clifford_unit(mu, A));
+
+       e = clifford_unit(nu, A) * clifford_unit(mu, A) * clifford_unit(nu_TOGGLE, A);
        
-       CHECK6(varidx,.toggle_variance())
+       result += check_equal_simplify_term(e,  2 * indexed(A_symm, sy_symm(), nu_TOGGLE, mu) *clifford_unit(nu, A)-A.trace()*clifford_unit(mu, A), mu);
 
-       return result;
-}
+       e = clifford_unit(nu, A) * clifford_unit(nu_TOGGLE, A)
+         * clifford_unit(mu, A) * clifford_unit(mu_TOGGLE, A);
+       result += check_equal_simplify(e, pow(A.trace(), 2) * dirac_ONE());
 
-static unsigned clifford_check6a(const matrix & A)
-{
-       matrix A_symm(4,4), A2(4, 4);
-       A_symm = A.add(A.transpose()).mul(half);
-       A2 = A_symm.mul(A_symm);
+       e = clifford_unit(mu, A) * clifford_unit(nu, A)
+         * clifford_unit(nu_TOGGLE, A) * clifford_unit(mu_TOGGLE, A);
+       result += check_equal_simplify(e, pow(A.trace(), 2)  * dirac_ONE());
+
+       e = clifford_unit(mu, A) * clifford_unit(nu, A)
+         * clifford_unit(mu_TOGGLE, A) * clifford_unit(nu_TOGGLE, A);
+
+       result += check_equal_simplify_term2(e, 2*indexed(A_symm, sy_symm(), nu_TOGGLE, mu_TOGGLE) * clifford_unit(mu, A) * clifford_unit(nu, A) - pow(A.trace(), 2)*dirac_ONE());
+
+       e = clifford_unit(mu_TOGGLE, A) * clifford_unit(nu, A)
+         * clifford_unit(mu, A) * clifford_unit(nu_TOGGLE, A);
+
+       result += check_equal_simplify_term2(e, 2*indexed(A_symm, nu, mu) * clifford_unit(mu_TOGGLE, A) * clifford_unit(nu_TOGGLE, A) - pow(A.trace(), 2)*dirac_ONE());
+
+       e = clifford_unit(nu_TOGGLE, A) * clifford_unit(rho_TOGGLE, A)
+         * clifford_unit(mu, A) * clifford_unit(rho, A) * clifford_unit(nu, A);
+       e = e.simplify_indexed().collect(clifford_unit(mu, A));
        
-       ex e, e1;
-       int result = 0;
+       result += check_equal_simplify_term(e, 4* indexed(A_symm, sy_symm(), nu_TOGGLE,  rho)*indexed(A_symm, sy_symm(), rho_TOGGLE, mu) *clifford_unit(nu, A) 
+                                           - 2*A.trace() * (clifford_unit(rho, A) * indexed(A_symm, sy_symm(), rho_TOGGLE, mu) 
+                                                            + clifford_unit(nu, A) * indexed(A_symm, sy_symm(), nu_TOGGLE, mu)) + pow(A.trace(),2)* clifford_unit(mu, A), mu);
 
-       CHECK6(idx,)
+       e = clifford_unit(nu_TOGGLE, A) * clifford_unit(rho, A)
+         * clifford_unit(mu, A) * clifford_unit(rho_TOGGLE, A) * clifford_unit(nu, A);
+       e = e.simplify_indexed().collect(clifford_unit(mu, A));
+       
+       result += check_equal_simplify_term(e, 4* indexed(A_symm, sy_symm(), nu_TOGGLE,  rho)*indexed(A_symm, sy_symm(), rho_TOGGLE, mu) *clifford_unit(nu, A) 
+                                           - 2*A.trace() * (clifford_unit(rho, A) * indexed(A_symm, sy_symm(), rho_TOGGLE, mu) 
+                                                            + clifford_unit(nu, A) * indexed(A_symm, sy_symm(), nu_TOGGLE, mu)) + pow(A.trace(),2)* clifford_unit(mu, A), mu);
+
+       e = clifford_unit(mu, A) * clifford_unit(nu, A) + clifford_unit(nu, A) * clifford_unit(mu, A);
+       result += check_equal(canonicalize_clifford(e), 2*dirac_ONE()*indexed(A_symm, sy_symm(), mu, nu));
+
+       e = (clifford_unit(mu, A) * clifford_unit(nu, A) * clifford_unit(lam, A)
+                + clifford_unit(nu, A) * clifford_unit(lam, A) * clifford_unit(mu, A)
+                + clifford_unit(lam, A) * clifford_unit(mu, A) * clifford_unit(nu, A)
+                - clifford_unit(nu, A) * clifford_unit(mu, A) * clifford_unit(lam, A)
+                - clifford_unit(lam, A) * clifford_unit(nu, A) * clifford_unit(mu, A)
+                - clifford_unit(mu, A) * clifford_unit(lam, A) * clifford_unit(nu, A)) / 6
+               + indexed(A_symm, sy_symm(), mu, nu) * clifford_unit(lam, A)
+               - indexed(A_symm, sy_symm(), mu, lam) * clifford_unit(nu, A)
+               + indexed(A_symm, sy_symm(), nu, lam) * clifford_unit(mu, A)
+               - clifford_unit(mu, A) * clifford_unit(nu, A) * clifford_unit(lam, A);
+       result += check_equal(canonicalize_clifford(e), 0);
 
+/* lst_to_clifford() and clifford_inverse()  check*/
+       realsymbol x("x"), y("y"), t("t"), z("z");
+
+       ex c = clifford_unit(nu, A, 1);
+       e = lst_to_clifford(lst(t, x, y, z), mu, A, 1) * lst_to_clifford(lst(1, 2, 3, 4), c);
+       e1 = clifford_inverse(e);
+       result += check_equal_simplify_term2((e*e1).simplify_indexed(), dirac_ONE(1));
+
+/* Moebius map (both forms) checks for symmetric metrics only */
+       matrix M1(2, 2),  M2(2, 2);
+       c = clifford_unit(nu, A);
+
+       e = clifford_moebius_map(0, dirac_ONE(), 
+                                                        dirac_ONE(), 0, lst(t, x, y, z), A); 
+/* this is just the inversion*/
+       M1 = 0, dirac_ONE(),
+               dirac_ONE(), 0;
+       e1 = clifford_moebius_map(M1, lst(t, x, y, z), A); 
+/* the inversion again*/
+       result += check_equal_lst(e, e1);
+
+       e1 = clifford_to_lst(clifford_inverse(lst_to_clifford(lst(t, x, y, z), mu, A)), c);
+       result += check_equal_lst(e, e1);
+
+       e = clifford_moebius_map(dirac_ONE(), lst_to_clifford(lst(1, 2, 3, 4), nu, A), 
+                                                        0, dirac_ONE(), lst(t, x, y, z), A); 
+/*this is just a shift*/
+       M2 = dirac_ONE(), lst_to_clifford(lst(1, 2, 3, 4), c),
+               0, dirac_ONE();
+       e1 = clifford_moebius_map(M2, lst(t, x, y, z), c); 
+/* the same shift*/
+       result += check_equal_lst(e, e1);
+
+       result += check_equal(e, lst(t+1, x+2, y+3, z+4));
+
+/* Check the group law for Moebius maps */
+       e = clifford_moebius_map(M1, ex_to<lst>(e1), c);
+/*composition of M1 and M2*/
+       e1 = clifford_moebius_map(M1.mul(M2), lst(t, x, y, z), c);
+/* the product M1*M2*/
+       result += check_equal_lst(e, e1);
        return result;
 }
 
@@ -545,45 +538,46 @@ unsigned exam_clifford()
        result += clifford_check5(); cout << '.' << flush;
 
        // anticommuting, symmetric examples
-       result += clifford_check6(ex_to<matrix>(diag_matrix(lst(-1, 1, 1, 1))))+clifford_check6a(ex_to<matrix>(diag_matrix(lst(-1, 1, 1, 1))));; cout << '.' << flush;
-       result += clifford_check6(ex_to<matrix>(diag_matrix(lst(-1, -1, -1, -1))))+clifford_check6a(ex_to<matrix>(diag_matrix(lst(-1, -1, -1, -1))));; cout << '.' << flush;
-       result += clifford_check6(ex_to<matrix>(diag_matrix(lst(-1, 1, 1, -1))))+clifford_check6a(ex_to<matrix>(diag_matrix(lst(-1, 1, 1, -1))));; cout << '.' << flush;
-       result += clifford_check6(ex_to<matrix>(diag_matrix(lst(-1, 0, 1, -1))))+clifford_check6a(ex_to<matrix>(diag_matrix(lst(-1, 0, 1, -1))));; cout << '.' << flush;
-       result += clifford_check6(ex_to<matrix>(diag_matrix(lst(-3, 0, 2, -1))))+clifford_check6a(ex_to<matrix>(diag_matrix(lst(-3, 0, 2, -1))));; cout << '.' << flush;
+       result += clifford_check6<varidx>(ex_to<matrix>(diag_matrix(lst(-1, 1, 1, 1))));
+       result += clifford_check6<idx>(ex_to<matrix>(diag_matrix(lst(-1, 1, 1, 1))));; cout << '.' << flush;
+       result += clifford_check6<varidx>(ex_to<matrix>(diag_matrix(lst(-1, -1, -1, -1))))+clifford_check6<idx>(ex_to<matrix>(diag_matrix(lst(-1, -1, -1, -1))));; cout << '.' << flush;
+       result += clifford_check6<idx>(ex_to<matrix>(diag_matrix(lst(-1, 1, 1, -1))))+clifford_check6<idx>(ex_to<matrix>(diag_matrix(lst(-1, 1, 1, -1))));; cout << '.' << flush;
+       result += clifford_check6<varidx>(ex_to<matrix>(diag_matrix(lst(-1, 0, 1, -1))))+clifford_check6<idx>(ex_to<matrix>(diag_matrix(lst(-1, 0, 1, -1))));; cout << '.' << flush;
+       result += clifford_check6<varidx>(ex_to<matrix>(diag_matrix(lst(-3, 0, 2, -1))))+clifford_check6<idx>(ex_to<matrix>(diag_matrix(lst(-3, 0, 2, -1))));; cout << '.' << flush;
 
        realsymbol s("s"), t("t"); // symbolic entries in matric
-       result += clifford_check6(ex_to<matrix>(diag_matrix(lst(-1, 1, s, t))))+clifford_check6a(ex_to<matrix>(diag_matrix(lst(-1, 1, s, t))));; cout << '.' << flush;
+       result += clifford_check6<varidx>(ex_to<matrix>(diag_matrix(lst(-1, 1, s, t))))+clifford_check6<idx>(ex_to<matrix>(diag_matrix(lst(-1, 1, s, t))));; cout << '.' << flush;
 
        matrix A(4, 4);
        A = 1, 0, 0, 0, // anticommuting, not symmetric, Tr=0
                0, -1, 0, 0,
                0, 0, 0, -1,
                0, 0, 1, 0; 
-       result += clifford_check6(A)+clifford_check6a(A);; cout << '.' << flush;
+       result += clifford_check6<varidx>(A)+clifford_check6<idx>(A);; cout << '.' << flush;
 
        A = 1, 0, 0, 0, // anticommuting, not symmetric, Tr=2
                0, 1, 0, 0,
                0, 0, 0, -1,
                0, 0, 1, 0; 
-       result += clifford_check6(A)+clifford_check6a(A);; cout << '.' << flush;
+       result += clifford_check6<varidx>(A)+clifford_check6<idx>(A);; cout << '.' << flush;
 
        A = 1, 0, 0, 0, // not anticommuting, symmetric, Tr=0
                0, -1, 0, 0,
                0, 0, 0, -1,
                0, 0, -1, 0; 
-       result += clifford_check6(A)+clifford_check6a(A);; cout << '.' << flush;
+       result += clifford_check6<varidx>(A)+clifford_check6<idx>(A);; cout << '.' << flush;
 
        A = 1, 0, 0, 0, // not anticommuting, symmetric, Tr=2
                0, 1, 0, 0,
                0, 0, 0, -1,
                0, 0, -1, 0; 
-       result += clifford_check6(A)+clifford_check6a(A);; cout << '.' << flush;
+       result += clifford_check6<varidx>(A)+clifford_check6<idx>(A);; cout << '.' << flush;
 
        A = 1, 1, 0, 0, // not anticommuting, not symmetric, Tr=4
                0, 1, 1, 0,
                0, 0, 1, 1,
                0, 0, 0, 1; 
-       result += clifford_check6(A)+clifford_check6a(A);; cout << '.' << flush;
+       result += clifford_check6<varidx>(A)+clifford_check6<idx>(A);; cout << '.' << flush;
 
        symbol dim("D");
        result += clifford_check7(minkmetric(), dim); cout << '.' << flush;

diff --git a/check/exam_paranoia.cpp b/check/exam_paranoia.cpp
index 7c82e48662d9181256a5ebc044edd2b501f3d3e4..0ab64910c437e7d8c51d773cdacc5637680307b5 100644 (file)
--- 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 d772b3f4743bd85ff8c127110f888c34801c0c6c..445c273ed62057b297843f880dbd1d3d78bd4363 100644 (file)
--- a/ginac/indexed.cpp
+++ b/ginac/indexed.cpp
@@ -639,60 +639,60 @@ bool reposition_dummy_indices(ex & e, exvector & variant_dummy_indices, exvector
 {
        bool something_changed = false;
 
-       // Canonicalize wrt indices that are dummies within e. I.e., their
-       // symbol occurs twice in an index of e. This is only done if there
-       // is a cyclic symmetry because in that case it may happen that after
-       // raising/lowering an index the indices get reshuffled by ::eval in
-       // such a way that the iterators no longer point to the right objects.
-       if (ex_to<symmetry>(ex_to<indexed>(e).get_symmetry()).has_cyclic()) {
-               // Get dummy pairs of varidxes within the indexed object in e.
-               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;
-                                               }
+       // 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;
                                }
+                               break;
                        }
                }
-               // 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;
+       }
+
+       // 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;
        }
+       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;
 
@@ -700,14 +700,20 @@ bool reposition_dummy_indices(ex & e, exvector & variant_dummy_indices, exvector
                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);
@@ -718,11 +724,8 @@ bool reposition_dummy_indices(ex & e, exvector & variant_dummy_indices, exvector
                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;
                        }
@@ -731,6 +734,9 @@ bool reposition_dummy_indices(ex & e, exvector & variant_dummy_indices, exvector
 next_index: ;
        }
 
+       if (something_changed)
+               e = ex_to<indexed>(e).thiscontainer(seq);
+
        return something_changed;
 }