Index renaming issues, sped up simplify_indexed, used defined NC-objects

[ginac.git] / ginac / indexed.cpp

diff --git a/ginac/indexed.cpp b/ginac/indexed.cpp
index ef904f7222bfc798420dee5cdb6bed6373426a26..64eba6d08cf1cccf6dbf455cbcce4c01356dd3e7 100644 (file)
--- a/ginac/indexed.cpp
+++ b/ginac/indexed.cpp
@@ -297,6 +297,9 @@ ex indexed::eval(int level) const
                return f * thiscontainer(v);
        }
 
                return f * thiscontainer(v);
        }
 

+       if(this->tinfo()==TINFO_indexed && seq.size()==1)
+               return base;
+

        // Canonicalize indices according to the symmetry properties
        if (seq.size() > 2) {
                exvector v = seq;
        // Canonicalize indices according to the symmetry properties
        if (seq.size() > 2) {
                exvector v = seq;


@@ -324,6 +327,14 @@ ex indexed::thiscontainer(std::auto_ptr<exvector> vp) const
        return indexed(ex_to<symmetry>(symtree), vp);
 }
 
        return indexed(ex_to<symmetry>(symtree), vp);
 }
 

+unsigned indexed::return_type() const
+{
+       if(is_a<matrix>(op(0)))
+               return return_types::commutative;
+       else
+               return op(0).return_type();
+}
+

 ex indexed::expand(unsigned options) const
 {
        GINAC_ASSERT(seq.size() > 0);
 ex indexed::expand(unsigned options) const
 {
        GINAC_ASSERT(seq.size() > 0);


@@ -532,6 +543,15 @@ exvector integral::get_free_indices() const
        return f.get_free_indices();
 }
 
        return f.get_free_indices();
 }
 

+template<class T> size_t number_of_type(const exvector&v)
+{
+       size_t number = 0;
+       for(exvector::const_iterator i=v.begin(); i!=v.end(); ++i)
+               if(is_exactly_a<T>(*i))
+                       ++number;
+       return number;
+}
+

 /** Rename dummy indices in an expression.
  *
  *  @param e Expression to work on
 /** Rename dummy indices in an expression.
  *
  *  @param e Expression to work on


@@ -540,10 +560,10 @@ exvector integral::get_free_indices() const
  *  @param global_dummy_indices The set of dummy indices that have appeared
  *    before and which we would like to use in "e", too. This gets updated
  *    by the function */
  *  @param global_dummy_indices The set of dummy indices that have appeared
  *    before and which we would like to use in "e", too. This gets updated
  *    by the function */

-static ex rename_dummy_indices(const ex & e, exvector & global_dummy_indices, exvector & local_dummy_indices)
+template<class T> static ex rename_dummy_indices(const ex & e, exvector & global_dummy_indices, exvector & local_dummy_indices)

 {
 {

-       size_t global_size = global_dummy_indices.size(),
-              local_size = local_dummy_indices.size();
+       size_t global_size = number_of_type<T>(global_dummy_indices),
+              local_size = number_of_type<T>(local_dummy_indices);

 
        // Any local dummy indices at all?
        if (local_size == 0)
 
        // Any local dummy indices at all?
        if (local_size == 0)


@@ -557,7 +577,7 @@ static ex rename_dummy_indices(const ex & e, exvector & global_dummy_indices, ex
                int remaining = local_size - global_size;
                exvector::const_iterator it = local_dummy_indices.begin(), itend = local_dummy_indices.end();
                while (it != itend && remaining > 0) {
                int remaining = local_size - global_size;
                exvector::const_iterator it = local_dummy_indices.begin(), itend = local_dummy_indices.end();
                while (it != itend && remaining > 0) {

-                       if (find_if(global_dummy_indices.begin(), global_dummy_indices.end(), bind2nd(op0_is_equal(), *it)) == global_dummy_indices.end()) {
+                       if (is_exactly_a<T>(*it) && find_if(global_dummy_indices.begin(), global_dummy_indices.end(), bind2nd(idx_is_equal_ignore_dim(), *it)) == global_dummy_indices.end()) {

                                global_dummy_indices.push_back(*it);
                                global_size++;
                                remaining--;
                                global_dummy_indices.push_back(*it);
                                global_size++;
                                remaining--;


@@ -575,11 +595,13 @@ static ex rename_dummy_indices(const ex & e, exvector & global_dummy_indices, ex
        exvector local_syms, global_syms;
        local_syms.reserve(local_size);
        global_syms.reserve(local_size);
        exvector local_syms, global_syms;
        local_syms.reserve(local_size);
        global_syms.reserve(local_size);

-       for (size_t i=0; i<local_size; i++)
-               local_syms.push_back(local_dummy_indices[i].op(0));
+       for (size_t i=0; local_syms.size()!=local_size; i++)
+               if(is_exactly_a<T>(local_dummy_indices[i]))
+                       local_syms.push_back(local_dummy_indices[i].op(0));

        shaker_sort(local_syms.begin(), local_syms.end(), ex_is_less(), ex_swap());
        shaker_sort(local_syms.begin(), local_syms.end(), ex_is_less(), ex_swap());

-       for (size_t i=0; i<local_size; i++) // don't use more global symbols than necessary
-               global_syms.push_back(global_dummy_indices[i].op(0));
+       for (size_t i=0; global_syms.size()!=local_size; i++) // don't use more global symbols than necessary
+               if(is_exactly_a<T>(global_dummy_indices[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
        shaker_sort(global_syms.begin(), global_syms.end(), ex_is_less(), ex_swap());
 
        // Remove common indices


@@ -704,6 +726,18 @@ static void product_to_exvector(const ex & e, exvector & v, bool & non_commutati
        }
 }
 
        }
 }
 

+template<class T> ex idx_symmetrization(const ex& r,const exvector& local_dummy_indices)
+{      exvector dummy_syms;
+       dummy_syms.reserve(r.nops());
+       for (exvector::const_iterator it = local_dummy_indices.begin(); it != local_dummy_indices.end(); ++it)
+                       if(is_exactly_a<T>(*it))
+                               dummy_syms.push_back(it->op(0));
+       if(dummy_syms.size() < 2)
+               return r;
+       ex q=symmetrize(r, dummy_syms);
+       return q;
+}
+

 /** 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)
 /** 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)


@@ -864,19 +898,26 @@ contraction_done:
        // The result should be symmetric with respect to exchange of dummy
        // indices, so if the symmetrization vanishes, the whole expression is
        // zero. This detects things like eps.i.j.k * p.j * p.k = 0.
        // The result should be symmetric with respect to exchange of dummy
        // indices, so if the symmetrization vanishes, the whole expression is
        // zero. This detects things like eps.i.j.k * p.j * p.k = 0.

-       if (local_dummy_indices.size() >= 2) {
-               exvector dummy_syms;
-               dummy_syms.reserve(local_dummy_indices.size());
-               for (exvector::const_iterator it = local_dummy_indices.begin(); it != local_dummy_indices.end(); ++it)
-                       dummy_syms.push_back(it->op(0));
-               if (symmetrize(r, dummy_syms).is_zero()) {
-                       free_indices.clear();
-                       return _ex0;
-               }
+       ex q = idx_symmetrization<idx>(r, local_dummy_indices);
+       if (q.is_zero()) {
+               free_indices.clear();
+               return _ex0;
+       }
+       q = idx_symmetrization<varidx>(q, local_dummy_indices);
+       if (q.is_zero()) {
+               free_indices.clear();
+               return _ex0;
+       }
+       q = idx_symmetrization<spinidx>(q, local_dummy_indices);
+       if (q.is_zero()) {
+               free_indices.clear();
+               return _ex0;

        }
 
        // Dummy index renaming
        }
 
        // Dummy index renaming

-       r = rename_dummy_indices(r, dummy_indices, local_dummy_indices);
+       r = rename_dummy_indices<idx>(r, dummy_indices, local_dummy_indices);
+       r = rename_dummy_indices<varidx>(r, dummy_indices, local_dummy_indices);
+       r = rename_dummy_indices<spinidx>(r, dummy_indices, local_dummy_indices);

 
        // Product of indexed object with a scalar?
        if (is_exactly_a<mul>(r) && r.nops() == 2
 
        // Product of indexed object with a scalar?
        if (is_exactly_a<mul>(r) && r.nops() == 2


@@ -943,6 +984,17 @@ public:
        }
 };
 
        }
 };
 

+bool hasindex(const ex &x, const ex &sym)
+{      
+       if(is_a<idx>(x) && x.op(0)==sym)
+               return true;
+       else
+               for(size_t i=0; i<x.nops(); ++i)
+                       if(hasindex(x.op(i), sym))
+                               return true;
+       return false;
+}
+

 /** Simplify indexed expression, return list of free indices. */
 ex simplify_indexed(const ex & e, exvector & free_indices, exvector & dummy_indices, const scalar_products & sp)
 {
 /** Simplify indexed expression, return list of free indices. */
 ex simplify_indexed(const ex & e, exvector & free_indices, exvector & dummy_indices, const scalar_products & sp)
 {


@@ -971,7 +1023,10 @@ ex simplify_indexed(const ex & e, exvector & free_indices, exvector & dummy_indi
                }
 
                // Rename the dummy indices
                }
 
                // Rename the dummy indices

-               return rename_dummy_indices(e_expanded, dummy_indices, local_dummy_indices);
+               e_expanded = rename_dummy_indices<idx>(e_expanded, dummy_indices, local_dummy_indices);
+               e_expanded = rename_dummy_indices<varidx>(e_expanded, dummy_indices, local_dummy_indices);
+               e_expanded = rename_dummy_indices<spinidx>(e_expanded, dummy_indices, local_dummy_indices);
+               return e_expanded;

        }
 
        // Simplification of sum = sum of simplifications, check consistency of
        }
 
        // Simplification of sum = sum of simplifications, check consistency of


@@ -1015,18 +1070,19 @@ ex simplify_indexed(const ex & e, exvector & free_indices, exvector & dummy_indi
                if (num_terms_orig < 2 || dummy_indices.size() < 2)
                        return sum;
 
                if (num_terms_orig < 2 || dummy_indices.size() < 2)
                        return sum;
 

-               // Yes, construct vector of all dummy index symbols
-               exvector dummy_syms;
-               dummy_syms.reserve(dummy_indices.size());
-               for (exvector::const_iterator it = dummy_indices.begin(); it != dummy_indices.end(); ++it)
-                       dummy_syms.push_back(it->op(0));
-

                // Chop the sum into terms and symmetrize each one over the dummy
                // indices
                std::vector<terminfo> terms;
                for (size_t i=0; i<sum.nops(); i++) {
                        const ex & term = sum.op(i);
                // Chop the sum into terms and symmetrize each one over the dummy
                // indices
                std::vector<terminfo> terms;
                for (size_t i=0; i<sum.nops(); i++) {
                        const ex & term = sum.op(i);

-                       ex term_symm = symmetrize(term, dummy_syms);
+                       exvector dummy_indices_of_term;
+                       dummy_indices_of_term.reserve(dummy_indices.size());
+                       for(exvector::iterator i=dummy_indices.begin(); i!=dummy_indices.end(); ++i)
+                               if(hasindex(term,i->op(0)))
+                                       dummy_indices_of_term.push_back(*i);
+                       ex term_symm = idx_symmetrization<idx>(term, dummy_indices_of_term);
+                       term_symm = idx_symmetrization<varidx>(term_symm, dummy_indices_of_term);
+                       term_symm = idx_symmetrization<spinidx>(term_symm, dummy_indices_of_term);

                        if (term_symm.is_zero())
                                continue;
                        terms.push_back(terminfo(term, term_symm));
                        if (term_symm.is_zero())
                                continue;
                        terms.push_back(terminfo(term, term_symm));


@@ -1318,9 +1374,9 @@ exvector get_all_dummy_indices(const ex & e)
        return v;
 }
 
        return v;
 }
 

-ex rename_dummy_indices_uniquely(const ex & a, const ex & b)
+ex rename_dummy_indices_uniquely(const exvector & va, const exvector & vb, const ex & b)

 {
 {

-       exvector va = get_all_dummy_indices(a), vb = get_all_dummy_indices(b), common_indices;
+       exvector common_indices;

        set_intersection(va.begin(), va.end(), vb.begin(), vb.end(), std::back_insert_iterator<exvector>(common_indices), ex_is_less());
        if (common_indices.empty()) {
                return b;
        set_intersection(va.begin(), va.end(), vb.begin(), vb.end(), std::back_insert_iterator<exvector>(common_indices), ex_is_less());
        if (common_indices.empty()) {
                return b;


@@ -1330,15 +1386,25 @@ ex rename_dummy_indices_uniquely(const ex & a, const ex & b)
                new_indices.reserve(2*common_indices.size());
                exvector::const_iterator ip = common_indices.begin(), ipend = common_indices.end();
                while (ip != ipend) {
                new_indices.reserve(2*common_indices.size());
                exvector::const_iterator ip = common_indices.begin(), ipend = common_indices.end();
                while (ip != ipend) {

-                       if (is_a<varidx>(*ip)) {
-                               varidx mu((new symbol)->setflag(status_flags::dynallocated), ex_to<varidx>(*ip).get_dim(), ex_to<varidx>(*ip).is_covariant());
-                               old_indices.push_back(*ip);
-                               new_indices.push_back(mu);
+                       ex newsym=(new symbol)->setflag(status_flags::dynallocated);
+                       ex newidx;
+                       if(is_exactly_a<spinidx>(*ip))
+                               newidx = (new spinidx(newsym, ex_to<spinidx>(*ip).get_dim(),
+                                               ex_to<spinidx>(*ip).is_covariant(),
+                                               ex_to<spinidx>(*ip).is_dotted()))
+                                       -> setflag(status_flags::dynallocated);
+                       else if (is_exactly_a<varidx>(*ip))
+                               newidx = (new varidx(newsym, ex_to<varidx>(*ip).get_dim(),
+                                               ex_to<varidx>(*ip).is_covariant()))
+                                       -> setflag(status_flags::dynallocated);
+                       else
+                               newidx = (new idx(newsym, ex_to<idx>(*ip).get_dim()))
+                                       -> setflag(status_flags::dynallocated);
+                       old_indices.push_back(*ip);
+                       new_indices.push_back(newidx);
+                       if(is_a<varidx>(*ip)) {

                                old_indices.push_back(ex_to<varidx>(*ip).toggle_variance());
                                old_indices.push_back(ex_to<varidx>(*ip).toggle_variance());

-                               new_indices.push_back(mu.toggle_variance());
-                       } else {
-                               old_indices.push_back(*ip);
-                               new_indices.push_back(idx((new symbol)->setflag(status_flags::dynallocated), ex_to<varidx>(*ip).get_dim()));
+                               new_indices.push_back(ex_to<varidx>(newidx).toggle_variance());

                        }
                        ++ip;
                }
                        }
                        ++ip;
                }


@@ -1346,6 +1412,15 @@ ex rename_dummy_indices_uniquely(const ex & a, const ex & b)
        }
 }
 
        }
 }
 

+ex rename_dummy_indices_uniquely(const ex & a, const ex & b)
+{
+       exvector va = get_all_dummy_indices(a);
+       exvector vb = get_all_dummy_indices(b);
+       sort(va.begin(), va.end(), ex_is_less());
+       sort(vb.begin(), vb.end(), ex_is_less());
+       return rename_dummy_indices_uniquely(va, vb, b);
+}
+

 ex expand_dummy_sum(const ex & e, bool subs_idx)
 {
        ex e_expanded = e.expand();
 ex expand_dummy_sum(const ex & e, bool subs_idx)
 {
        ex e_expanded = e.expand();