Synced to HEAD

[ginac.git] / ginac / indexed.cpp

diff --git a/ginac/indexed.cpp b/ginac/indexed.cpp
index 0bb34d8a578a4f42315aaf7b0d13a854b74cefe8..d693373cae62c6cd622495d05722aa9ce7065d99 100644 (file)
--- a/ginac/indexed.cpp
+++ b/ginac/indexed.cpp
@@ -48,7 +48,7 @@ GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(indexed, exprseq,
 // default constructor
 //////////
 
-indexed::indexed() : symtree(sy_none())
+indexed::indexed() : symtree(not_symmetric())
 {
        tinfo_key = TINFO_indexed;
 }
@@ -57,31 +57,31 @@ indexed::indexed() : symtree(sy_none())
 // other constructors
 //////////
 
-indexed::indexed(const ex & b) : inherited(b), symtree(sy_none())
+indexed::indexed(const ex & b) : inherited(b), symtree(not_symmetric())
 {
        tinfo_key = TINFO_indexed;
        validate();
 }
 
-indexed::indexed(const ex & b, const ex & i1) : inherited(b, i1), symtree(sy_none())
+indexed::indexed(const ex & b, const ex & i1) : inherited(b, i1), symtree(not_symmetric())
 {
        tinfo_key = TINFO_indexed;
        validate();
 }
 
-indexed::indexed(const ex & b, const ex & i1, const ex & i2) : inherited(b, i1, i2), symtree(sy_none())
+indexed::indexed(const ex & b, const ex & i1, const ex & i2) : inherited(b, i1, i2), symtree(not_symmetric())
 {
        tinfo_key = TINFO_indexed;
        validate();
 }
 
-indexed::indexed(const ex & b, const ex & i1, const ex & i2, const ex & i3) : inherited(b, i1, i2, i3), symtree(sy_none())
+indexed::indexed(const ex & b, const ex & i1, const ex & i2, const ex & i3) : inherited(b, i1, i2, i3), symtree(not_symmetric())
 {
        tinfo_key = TINFO_indexed;
        validate();
 }
 
-indexed::indexed(const ex & b, const ex & i1, const ex & i2, const ex & i3, const ex & i4) : inherited(b, i1, i2, i3, i4), symtree(sy_none())
+indexed::indexed(const ex & b, const ex & i1, const ex & i2, const ex & i3, const ex & i4) : inherited(b, i1, i2, i3, i4), symtree(not_symmetric())
 {
        tinfo_key = TINFO_indexed;
        validate();
@@ -105,7 +105,7 @@ indexed::indexed(const ex & b, const symmetry & symm, const ex & i1, const ex &
        validate();
 }
 
-indexed::indexed(const ex & b, const exvector & v) : inherited(b), symtree(sy_none())
+indexed::indexed(const ex & b, const exvector & v) : inherited(b), symtree(not_symmetric())
 {
        seq.insert(seq.end(), v.begin(), v.end());
        tinfo_key = TINFO_indexed;
@@ -152,7 +152,7 @@ indexed::indexed(const archive_node &n, lst &sym_lst) : inherited(n, sym_lst)
                                symtree = sy_anti();
                                break;
                        default:
-                               symtree = sy_none();
+                               symtree = not_symmetric();
                                break;
                }
                const_cast<symmetry &>(ex_to<symmetry>(symtree)).validate(seq.size() - 1);
@@ -325,20 +325,24 @@ ex indexed::expand(unsigned options) const
 {
        GINAC_ASSERT(seq.size() > 0);
 
-       if ((options & expand_options::expand_indexed) && is_exactly_a<add>(seq[0])) {
-
-               // expand_indexed expands (a+b).i -> a.i + b.i
-               const ex & base = seq[0];
-               ex sum = _ex0;
-               for (size_t i=0; i<base.nops(); i++) {
+       if (options & expand_options::expand_indexed) {
+               ex newbase = seq[0].expand(options);
+               if (is_exactly_a<add>(newbase)) {
+                       ex sum = _ex0;
+                       for (size_t i=0; i<newbase.nops(); i++) {
+                               exvector s = seq;
+                               s[0] = newbase.op(i);
+                               sum += thiscontainer(s).expand(options);
+                       }
+                       return sum;
+               }
+               if (!are_ex_trivially_equal(newbase, seq[0])) {
                        exvector s = seq;
-                       s[0] = base.op(i);
-                       sum += thiscontainer(s).expand();
+                       s[0] = newbase;
+                       return ex_to<indexed>(thiscontainer(s)).inherited::expand(options);
                }
-               return sum;
-
-       } else
-               return inherited::expand(options);
+       }
+       return inherited::expand(options);
 }
 
 //////////
@@ -495,10 +499,27 @@ exvector ncmul::get_free_indices() const
        return free_indices;
 }
 
+struct is_summation_idx : public std::unary_function<ex, bool> {
+       bool operator()(const ex & e)
+       {
+               return is_dummy_pair(e, e);
+       }
+};
+
 exvector power::get_free_indices() const
 {
-       // Return free indices of basis
-       return basis.get_free_indices();
+       // Get free indices of basis
+       exvector basis_indices = basis.get_free_indices();
+
+       if (exponent.info(info_flags::even)) {
+               // If the exponent is an even number, then any "free" index that
+               // forms a dummy pair with itself is actually a summation index
+               exvector really_free;
+               std::remove_copy_if(basis_indices.begin(), basis_indices.end(),
+                                   std::back_inserter(really_free), is_summation_idx());
+               return really_free;
+       } else
+               return basis_indices;
 }
 
 /** Rename dummy indices in an expression.
@@ -1107,6 +1128,7 @@ ex simplify_indexed(const ex & e, exvector & free_indices, exvector & dummy_indi
  *  performs contraction of dummy indices where possible and checks whether
  *  the free indices in sums are consistent.
  *
+ *  @param options Simplification options (currently unused)
  *  @return simplified expression */
 ex ex::simplify_indexed(unsigned options) const
 {
@@ -1121,6 +1143,7 @@ ex ex::simplify_indexed(unsigned options) const
  *  scalar products by known values if desired.
  *
  *  @param sp Scalar products to be replaced automatically
+ *  @param options Simplification options (currently unused)
  *  @return simplified expression */
 ex ex::simplify_indexed(const scalar_products & sp, unsigned options) const
 {