Improved clifford_moebius_map() [V.Kisil]

[ginac.git] / ginac / indexed.cpp

diff --git a/ginac/indexed.cpp b/ginac/indexed.cpp
index bd41f44b0b38a271d243564c945d674bea9acaa3..b3c18ba51989ba37bdebb03b111f007858721ec1 100644 (file)
--- a/ginac/indexed.cpp
+++ b/ginac/indexed.cpp
@@ -3,7 +3,7 @@
  *  Implementation of GiNaC's indexed expressions. */
 
 /*
- *  GiNaC Copyright (C) 1999-2006 Johannes Gutenberg University Mainz, Germany
+ *  GiNaC Copyright (C) 1999-2007 Johannes Gutenberg University Mainz, Germany
  *
  *  This program is free software; you can redistribute it and/or modify
  *  it under the terms of the GNU General Public License as published by
@@ -23,6 +23,7 @@
 #include <iostream>
 #include <sstream>
 #include <stdexcept>
+#include <limits>
 
 #include "indexed.h"
 #include "idx.h"
@@ -306,7 +307,7 @@ ex indexed::eval(int level) const
                exvector v = seq;
                GINAC_ASSERT(is_exactly_a<symmetry>(symtree));
                int sig = canonicalize(v.begin() + 1, ex_to<symmetry>(symtree));
-               if (sig != INT_MAX) {
+               if (sig != std::numeric_limits<int>::max()) {
                        // Something has changed while sorting indices, more evaluations later
                        if (sig == 0)
                                return _ex0;
@@ -639,10 +640,60 @@ bool reposition_dummy_indices(ex & e, exvector & variant_dummy_indices, exvector
 {
        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;
 
@@ -650,14 +701,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);
@@ -668,11 +725,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;
                        }
@@ -681,6 +735,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;
 }