]> www.ginac.de Git - ginac.git/blobdiff - ginac/indexed.cpp
simplify_indexed() raises/lowers dummy indices to canonicalize their variance
[ginac.git] / ginac / indexed.cpp
index a89fb1011ed42eb2d5882e8e79ba4c101f3cbb55..eec7b11fe6fa4cd8057d1059240172b7de31f590 100644 (file)
@@ -3,7 +3,7 @@
  *  Implementation of GiNaC's indexed expressions. */
 
 /*
- *  GiNaC Copyright (C) 1999-2001 Johannes Gutenberg University Mainz, Germany
+ *  GiNaC Copyright (C) 1999-2002 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
@@ -20,8 +20,8 @@
  *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
  */
 
+#include <iostream>
 #include <stdexcept>
-#include <algorithm>
 
 #include "indexed.h"
 #include "idx.h"
@@ -29,6 +29,7 @@
 #include "mul.h"
 #include "ncmul.h"
 #include "power.h"
+#include "relational.h"
 #include "symmetry.h"
 #include "lst.h"
 #include "print.h"
@@ -274,24 +275,6 @@ ex indexed::eval(int level) const
        return ex_to<basic>(base).eval_indexed(*this);
 }
 
-int indexed::degree(const ex & s) const
-{
-       return is_equal(ex_to<basic>(s)) ? 1 : 0;
-}
-
-int indexed::ldegree(const ex & s) const
-{
-       return is_equal(ex_to<basic>(s)) ? 1 : 0;
-}
-
-ex indexed::coeff(const ex & s, int n) const
-{
-       if (is_equal(ex_to<basic>(s)))
-               return n==1 ? _ex1 : _ex0;
-       else
-               return n==0 ? ex(*this) : _ex0;
-}
-
 ex indexed::thisexprseq(const exvector & v) const
 {
        return indexed(ex_to<symmetry>(symtree), v);
@@ -346,9 +329,10 @@ void indexed::printindices(const print_context & c, unsigned level) const
 
                        while (it != itend) {
                                bool cur_covariant = (is_ex_of_type(*it, varidx) ? ex_to<varidx>(*it).is_covariant() : true);
-                               if (first || cur_covariant != covariant) {
+                               if (first || cur_covariant != covariant) { // Variance changed
+                                       // The empty {} prevents indices from ending up on top of each other
                                        if (!first)
-                                               c.s << "}";
+                                               c.s << "}{}";
                                        covariant = cur_covariant;
                                        if (covariant)
                                                c.s << "_{";
@@ -539,7 +523,6 @@ static ex rename_dummy_indices(const ex & e, exvector & global_dummy_indices, ex
                        }
                        it++;
                }
-               shaker_sort(global_dummy_indices.begin(), global_dummy_indices.end(), ex_is_less(), ex_swap());
 
                // If this is the first set of local indices, do nothing
                if (old_global_size == 0)
@@ -554,6 +537,7 @@ static ex rename_dummy_indices(const ex & e, exvector & global_dummy_indices, ex
        shaker_sort(local_syms.begin(), local_syms.end(), ex_is_less(), ex_swap());
        for (unsigned i=0; i<global_size; i++)
                global_syms.push_back(global_dummy_indices[i].op(0));
+       shaker_sort(global_syms.begin(), global_syms.end(), ex_is_less(), ex_swap());
 
        // Remove common indices
        exlist local_uniq, global_uniq;
@@ -570,6 +554,14 @@ static ex rename_dummy_indices(const ex & e, exvector & global_dummy_indices, ex
        }
 }
 
+/* Ordering that only compares the base expressions of indexed objects. */
+struct ex_base_is_less : public std::binary_function<ex, ex, bool> {
+       bool operator() (const ex &lh, const ex &rh) const
+       {
+               return (is_a<indexed>(lh) ? lh.op(0) : lh).compare(is_a<indexed>(rh) ? rh.op(0) : rh) < 0;
+       }
+};
+
 /** Simplify product of indexed expressions (commutative, noncommutative and
  *  simple squares), return list of free indices. */
 ex simplify_indexed_product(const ex & e, exvector & free_indices, exvector & dummy_indices, const scalar_products & sp)
@@ -687,8 +679,7 @@ contraction_done:
        // Find free indices (concatenate them all and call find_free_and_dummy())
        // and all dummy indices that appear
        exvector un, individual_dummy_indices;
-       it1 = v.begin(); itend = v.end();
-       while (it1 != itend) {
+       for (it1 = v.begin(), itend = v.end(); it1 != itend; ++it1) {
                exvector free_indices_of_factor;
                if (is_ex_of_type(*it1, indexed)) {
                        exvector dummy_indices_of_factor;
@@ -697,12 +688,76 @@ contraction_done:
                } else
                        free_indices_of_factor = it1->get_free_indices();
                un.insert(un.end(), free_indices_of_factor.begin(), free_indices_of_factor.end());
-               it1++;
        }
        exvector local_dummy_indices;
        find_free_and_dummy(un, free_indices, local_dummy_indices);
        local_dummy_indices.insert(local_dummy_indices.end(), individual_dummy_indices.begin(), individual_dummy_indices.end());
 
+       // Filter out the dummy indices with variance
+       exvector variant_dummy_indices;
+       for (it1 = local_dummy_indices.begin(), itend = local_dummy_indices.end(); it1 != itend; ++it1) {
+               if (is_exactly_a<varidx>(*it1))
+                       variant_dummy_indices.push_back(*it1);
+       }
+
+       // Any indices with variance present at all?
+       if (!variant_dummy_indices.empty()) {
+
+               // Yes, bring the product into a canonical order that only depends on
+               // the base expressions of indexed objects
+               if (!non_commutative)
+                       std::sort(v.begin(), v.end(), ex_base_is_less());
+
+               exvector moved_indices;
+
+               // Iterate over all indexed objects in the product
+               for (it1 = v.begin(), itend = v.end(); it1 != itend; ++it1) {
+                       if (!is_ex_of_type(*it1, indexed))
+                               continue;
+
+                       ex new_it1;
+                       bool it1_dirty = false; // It this is true, then new_it1 holds a new value for *it1
+
+                       // If a dummy index is encountered for the first time in the
+                       // product, pull it up, otherwise, pull it down
+                       exvector::iterator it2, it2end;
+                       for (it2 = const_cast<indexed &>(ex_to<indexed>(*it1)).seq.begin(), it2end = const_cast<indexed &>(ex_to<indexed>(*it1)).seq.end(); 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()) {
+                                                       new_it1 = (it1_dirty ? new_it1 : *it1).subs(*it2 == ex_to<varidx>(*it2).toggle_variance());
+                                                       it1_dirty = true;
+                                                       something_changed = true;
+                                               }
+                                               moved_indices.push_back(*vit);
+                                               variant_dummy_indices.erase(vit);
+                                               goto next_index;
+                                       }
+                               }
+
+                               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()) {
+                                                       new_it1 = (it1_dirty ? new_it1 : *it1).subs(*it2 == ex_to<varidx>(*it2).toggle_variance());
+                                                       it1_dirty = true;
+                                                       something_changed = true;
+                                               }
+                                               goto next_index;
+                                       }
+                               }
+
+next_index:            ;
+                       }
+
+                       if (it1_dirty)
+                               *it1 = new_it1;
+               }
+       }
+
        ex r;
        if (something_changed)
                r = non_commutative ? ex(ncmul(v, true)) : ex(mul(v));