+
+ // 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;
+ ex term = simplify_indexed(e_expanded.op(i), free_indices_of_term, dummy_indices, 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 = ex_to<basic>(sum.op(0)).add_indexed(sum, term);
+ else
+ sum += term;
+ }
+ }
+ }
+
+ if (sum.is_zero()) {
+ free_indices.clear();
+ return sum;
+ }
+
+ // Symmetrizing over the dummy indices may cancel terms
+ int num_terms_orig = (is_a<add>(sum) ? sum.nops() : 1);
+ if (num_terms_orig > 1 && dummy_indices.size() >= 2) {
+ lst dummy_syms;
+ for (int i=0; i<dummy_indices.size(); i++)
+ dummy_syms.append(dummy_indices[i].op(0));
+ ex sum_symm = sum.symmetrize(dummy_syms);
+ if (sum_symm.is_zero()) {
+ free_indices.clear();
+ return _ex0;
+ }
+ int num_terms = (is_a<add>(sum_symm) ? sum_symm.nops() : 1);
+ if (num_terms < num_terms_orig)
+ return sum_symm;
+ }
+
+ return sum;
+ }
+
+ // Simplification of products
+ if (is_ex_exactly_of_type(e_expanded, mul)
+ || is_ex_exactly_of_type(e_expanded, ncmul)
+ || (is_ex_exactly_of_type(e_expanded, power) && is_ex_of_type(e_expanded.op(0), indexed) && e_expanded.op(1).is_equal(_ex2)))
+ return simplify_indexed_product(e_expanded, free_indices, dummy_indices, sp);
+
+ // Cannot do anything
+ free_indices.clear();
+ return e_expanded;
+}
+
+/** Simplify/canonicalize expression containing indexed objects. This
+ * performs contraction of dummy indices where possible and checks whether
+ * the free indices in sums are consistent.
+ *
+ * @return simplified expression */
+ex ex::simplify_indexed(void) const
+{
+ exvector free_indices, dummy_indices;
+ scalar_products sp;
+ return GiNaC::simplify_indexed(*this, free_indices, dummy_indices, sp);
+}
+
+/** Simplify/canonicalize expression containing indexed objects. This
+ * performs contraction of dummy indices where possible, checks whether
+ * the free indices in sums are consistent, and automatically replaces
+ * scalar products by known values if desired.
+ *
+ * @param sp Scalar products to be replaced automatically
+ * @return simplified expression */
+ex ex::simplify_indexed(const scalar_products & sp) const
+{
+ exvector free_indices, dummy_indices;
+ return GiNaC::simplify_indexed(*this, free_indices, dummy_indices, sp);
+}
+
+/** Symmetrize expression over its free indices. */
+ex ex::symmetrize(void) const
+{
+ return GiNaC::symmetrize(*this, get_free_indices());
+}
+
+/** Antisymmetrize expression over its free indices. */
+ex ex::antisymmetrize(void) const
+{
+ return GiNaC::antisymmetrize(*this, get_free_indices());
+}
+
+/** Symmetrize expression by cyclic permutation over its free indices. */
+ex ex::symmetrize_cyclic(void) const
+{
+ return GiNaC::symmetrize_cyclic(*this, get_free_indices());
+}
+
+//////////
+// helper classes
+//////////
+
+void scalar_products::add(const ex & v1, const ex & v2, const ex & sp)
+{
+ spm[make_key(v1, v2)] = sp;
+}
+
+void scalar_products::add_vectors(const lst & l)
+{
+ // Add all possible pairs of products
+ unsigned num = l.nops();
+ for (unsigned i=0; i<num; i++) {
+ ex a = l.op(i);
+ for (unsigned j=0; j<num; j++) {
+ ex b = l.op(j);
+ add(a, b, a*b);
+ }
+ }
+}
+
+void scalar_products::clear(void)
+{
+ spm.clear();
+}
+
+/** Check whether scalar product pair is defined. */
+bool scalar_products::is_defined(const ex & v1, const ex & v2) const
+{
+ return spm.find(make_key(v1, v2)) != spm.end();
+}
+
+/** Return value of defined scalar product pair. */
+ex scalar_products::evaluate(const ex & v1, const ex & v2) const
+{
+ return spm.find(make_key(v1, v2))->second;
+}
+
+void scalar_products::debugprint(void) const
+{
+ std::cerr << "map size=" << spm.size() << std::endl;
+ spmap::const_iterator i = spm.begin(), end = spm.end();
+ while (i != end) {
+ const spmapkey & k = i->first;
+ std::cerr << "item key=(" << k.first << "," << k.second;
+ std::cerr << "), value=" << i->second << std::endl;
+ ++i;
+ }
+}
+
+/** Make key from object pair. */
+spmapkey scalar_products::make_key(const ex & v1, const ex & v2)
+{
+ // If indexed, extract base objects
+ ex s1 = is_ex_of_type(v1, indexed) ? v1.op(0) : v1;
+ ex s2 = is_ex_of_type(v2, indexed) ? v2.op(0) : v2;
+
+ // Enforce canonical order in pair
+ if (s1.compare(s2) > 0)
+ return spmapkey(s2, s1);
+ else
+ return spmapkey(s1, s2);