- subs() can be used to substitute functions, tensors and indexed objects

[ginac.git] / ginac / idx.cpp

diff --git a/ginac/idx.cpp b/ginac/idx.cpp
index 89373cc465ad4a770607cc98de49f55d95c563b5..0a8686de4593a06e61b2b252226e0809026c5edb 100644 (file)
--- a/ginac/idx.cpp
+++ b/ginac/idx.cpp
@@ -29,6 +29,8 @@
 #include "utils.h"
 #include "debugmsg.h"
 
+#include "exprseq.h" // !!
+
 namespace GiNaC {
 
 GINAC_IMPLEMENT_REGISTERED_CLASS(idx, basic)
@@ -62,17 +64,8 @@ void varidx::copy(const varidx & other)
        covariant = other.covariant;
 }
 
-void idx::destroy(bool call_parent)
-{
-       if (call_parent)
-               inherited::destroy(call_parent);
-}
-
-void varidx::destroy(bool call_parent)
-{
-       if (call_parent)
-               inherited::destroy(call_parent);
-}
+DEFAULT_DESTROY(idx)
+DEFAULT_DESTROY(varidx)
 
 //////////
 // other constructors
@@ -109,16 +102,6 @@ varidx::varidx(const archive_node &n, const lst &sym_lst) : inherited(n, sym_lst
        n.find_bool("covariant", covariant);
 }
 
-ex idx::unarchive(const archive_node &n, const lst &sym_lst)
-{
-       return (new idx(n, sym_lst))->setflag(status_flags::dynallocated);
-}
-
-ex varidx::unarchive(const archive_node &n, const lst &sym_lst)
-{
-       return (new varidx(n, sym_lst))->setflag(status_flags::dynallocated);
-}
-
 void idx::archive(archive_node &n) const
 {
        inherited::archive(n);
@@ -132,6 +115,9 @@ void varidx::archive(archive_node &n) const
        n.add_bool("covariant", covariant);
 }
 
+DEFAULT_UNARCHIVE(idx)
+DEFAULT_UNARCHIVE(varidx)
+
 //////////
 // functions overriding virtual functions from bases classes
 //////////
@@ -198,6 +184,18 @@ bool idx::info(unsigned inf) const
        return inherited::info(inf);
 }
 
+unsigned idx::nops() const
+{
+       // don't count the dimension as that is not really a sub-expression
+       return 1;
+}
+
+ex & idx::let_op(int i)
+{
+       GINAC_ASSERT(i == 0);
+       return value;
+}
+
 /** Returns order relation between two indices of the same type. The order
  *  must be such that dummy indices lie next to each other. */
 int idx::compare_same_type(const basic & other) const
@@ -241,6 +239,7 @@ ex idx::subs(const lst & ls, const lst & lr) const
                        // Otherwise substitute value
                        idx *i_copy = static_cast<idx *>(duplicate());
                        i_copy->value = lr.op(i);
+                       i_copy->clearflag(status_flags::hash_calculated);
                        return i_copy->setflag(status_flags::dynallocated);
                }
        }
@@ -252,6 +251,7 @@ ex idx::subs(const lst & ls, const lst & lr) const
 
        idx *i_copy = static_cast<idx *>(duplicate());
        i_copy->value = subsed_value;
+       i_copy->clearflag(status_flags::hash_calculated);
        return i_copy->setflag(status_flags::dynallocated);
 }
 
@@ -321,4 +321,90 @@ bool is_dummy_pair(const ex & e1, const ex & e2)
        return is_dummy_pair(ex_to_idx(e1), ex_to_idx(e2));
 }
 
+/** Bring a vector of indices into a canonic order. Dummy indices will lie
+ *  next to each other after the sorting. */
+static void sort_index_vector(exvector &v)
+{
+       // Nothing to sort if less than 2 elements
+       if (v.size() < 2)
+               return;
+
+       // Simple bubble sort algorithm should be sufficient for the small
+       // number of indices expected
+       exvector::iterator it1 = v.begin(), itend = v.end(), next_to_last_idx = itend - 1;
+       while (it1 != next_to_last_idx) {
+               exvector::iterator it2 = it1 + 1;
+               while (it2 != itend) {
+                       if (it1->compare(*it2) > 0)
+                               it1->swap(*it2);
+                       it2++;
+               }
+               it1++;
+       }
+}
+
+
+void find_free_and_dummy(exvector::const_iterator it, exvector::const_iterator itend, exvector & out_free, exvector & out_dummy)
+{
+       out_free.clear();
+       out_dummy.clear();
+
+       // No indices? Then do nothing
+       if (it == itend)
+               return;
+
+       // Only one index? Then it is a free one if it's not numeric
+       if (itend - it == 1) {
+               if (ex_to_idx(*it).is_symbolic())
+                       out_free.push_back(*it);
+               return;
+       }
+
+       // Sort index vector. This will cause dummy indices come to lie next
+       // to each other (because the sort order is defined to guarantee this).
+       exvector v(it, itend);
+       sort_index_vector(v);
+
+       // Find dummy pairs and free indices
+       it = v.begin(); itend = v.end();
+       exvector::const_iterator last = it++;
+       while (it != itend) {
+               if (is_dummy_pair(*it, *last)) {
+                       out_dummy.push_back(*last);
+                       it++;
+                       if (it == itend)
+                               return;
+               } else {
+                       if (!it->is_equal(*last) && ex_to_idx(*last).is_symbolic())
+                               out_free.push_back(*last);
+               }
+               last = it++;
+       }
+       if (ex_to_idx(*last).is_symbolic())
+               out_free.push_back(*last);
+}
+
+exvector index_set_difference(const exvector & set1, const exvector & set2)
+{
+       exvector ret;
+
+       exvector::const_iterator ait = set1.begin(), aitend = set1.end();
+       while (ait != aitend) {
+               exvector::const_iterator bit = set2.begin(), bitend = set2.end();
+               bool found = false;
+               while (bit != bitend) {
+                       if (ait->is_equal(*bit)) {
+                               found = true;
+                               break;
+                       }
+                       bit++;
+               }
+               if (!found)
+                       ret.push_back(*ait);
+               ait++;
+       }
+
+       return ret;
+}
+
 } // namespace GiNaC