if (level > 1)
return indexed(symmetry, evalchildren(level));
+ const ex &base = seq[0];
+
// If the base object is 0, the whole object is 0
- if (seq[0].is_zero())
+ if (base.is_zero())
return _ex0();
+ // If the base object is a product, pull out the numeric factor
+ if (is_ex_exactly_of_type(base, mul) && is_ex_exactly_of_type(base.op(base.nops() - 1), numeric)) {
+ exvector v = seq;
+ ex f = ex_to_numeric(base.op(base.nops() - 1));
+ v[0] = seq[0] / f;
+ return f * thisexprseq(v);
+ }
+
// Canonicalize indices according to the symmetry properties
if (seq.size() > 2 && (symmetry != unknown && symmetry != mixed)) {
exvector v = seq;
}
// Let the class of the base object perform additional evaluations
- return seq[0].bp->eval_indexed(*this);
+ return base.bp->eval_indexed(*this);
}
ex indexed::thisexprseq(const exvector & v) const
}
}
+ // Contraction of symmetric with antisymmetric object is zero
+ if ((ex_to_indexed(*it1).symmetry == indexed::symmetric &&
+ ex_to_indexed(*it2).symmetry == indexed::antisymmetric
+ || ex_to_indexed(*it1).symmetry == indexed::antisymmetric &&
+ ex_to_indexed(*it2).symmetry == indexed::symmetric)
+ && dummy.size() > 1) {
+ free_indices.clear();
+ return _ex0();
+ }
+
// Try to contract the first one with the second one
bool contracted = it1->op(0).bp->contract_with(it1, it2, v);
if (!contracted) {
}
find_free_and_dummy(un.begin(), un.end(), free_indices, dummy_indices);
+ ex r;
if (something_changed) {
if (non_commutative)
- return ncmul(v);
+ r = ncmul(v);
else
- return mul(v);
+ r = mul(v);
} else
- return e;
+ r = e;
+
+ // Product of indexed object with a scalar?
+ if (is_ex_exactly_of_type(r, mul) && r.nops() == 2
+ && is_ex_exactly_of_type(r.op(1), numeric) && is_ex_of_type(r.op(0), indexed))
+ return r.op(0).op(0).bp->scalar_mul_indexed(r.op(0), ex_to_numeric(r.op(1)));
+ else
+ return r;
}
/** Simplify indexed expression, return list of free indices. */
// Simplification of sum = sum of simplifications, check consistency of
// free indices in each term
if (is_ex_exactly_of_type(e_expanded, add)) {
+ bool first = true;
ex sum = _ex0();
+ free_indices.clear();
for (unsigned i=0; i<e_expanded.nops(); i++) {
exvector free_indices_of_term;
- sum += simplify_indexed(e_expanded.op(i), free_indices_of_term, sp);
- if (i == 0)
- free_indices = free_indices_of_term;
- else if (!indices_consistent(free_indices, free_indices_of_term))
- throw (std::runtime_error("simplify_indexed: inconsistent indices in sum"));
+ ex term = simplify_indexed(e_expanded.op(i), free_indices_of_term, sp);
+ if (!term.is_zero()) {
+ if (first) {
+ free_indices = free_indices_of_term;
+ sum = term;
+ first = false;
+ } else {
+ if (!indices_consistent(free_indices, free_indices_of_term))
+ throw (std::runtime_error("simplify_indexed: inconsistent indices in sum"));
+ if (is_ex_of_type(sum, indexed) && is_ex_of_type(term, indexed))
+ sum = sum.op(0).bp->add_indexed(sum, term);
+ else
+ sum += term;
+ }
+ }
}
return sum;