X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=blobdiff_plain;f=ginac%2Findexed.cpp;h=6fbe40b97cb726eaad9a24d395066d0511d6992a;hp=4fa4256abfb286b86c95d427efa9a5c8de27c4c8;hb=e455b8fcbb2788d437a7d212811bb5ba171e0cd4;hpb=fd097c33e414f72c77376ae7ad4fbd630b3ccce6

diff --git a/ginac/indexed.cpp b/ginac/indexed.cpp
index 4fa4256a..6fbe40b9 100644
--- a/ginac/indexed.cpp
+++ b/ginac/indexed.cpp
@@ -28,6 +28,8 @@
 #include "mul.h"
 #include "ncmul.h"
 #include "power.h"
+#include "lst.h"
+#include "print.h"
 #include "archive.h"
 #include "utils.h"
 #include "debugmsg.h"
@@ -52,11 +54,7 @@ void indexed::copy(const indexed & other)
 	symmetry = other.symmetry;
 }
 
-void indexed::destroy(bool call_parent)
-{
-	if (call_parent)
-		inherited::destroy(call_parent);
-}
+DEFAULT_DESTROY(indexed)
 
 //////////
 // other constructors
@@ -159,7 +157,6 @@ indexed::indexed(symmetry_type symm, exvector * vp) : inherited(vp), symmetry(sy
 // archiving
 //////////
 
-/** Construct object from archive_node. */
 indexed::indexed(const archive_node &n, const lst &sym_lst) : inherited(n, sym_lst)
 {
 	debugmsg("indexed constructor from archive_node", LOGLEVEL_CONSTRUCT);
@@ -168,59 +165,50 @@ indexed::indexed(const archive_node &n, const lst &sym_lst) : inherited(n, sym_l
 		throw (std::runtime_error("unknown indexed symmetry type in archive"));
 }
 
-/** Unarchive the object. */
-ex indexed::unarchive(const archive_node &n, const lst &sym_lst)
-{
-	return (new indexed(n, sym_lst))->setflag(status_flags::dynallocated);
-}
-
-/** Archive the object. */
 void indexed::archive(archive_node &n) const
 {
 	inherited::archive(n);
 	n.add_unsigned("symmetry", symmetry);
 }
 
+DEFAULT_UNARCHIVE(indexed)
+
 //////////
 // functions overriding virtual functions from bases classes
 //////////
 
-void indexed::printraw(std::ostream & os) const
+void indexed::print(const print_context & c, unsigned level) const
 {
-	debugmsg("indexed printraw", LOGLEVEL_PRINT);
+	debugmsg("indexed print", LOGLEVEL_PRINT);
 	GINAC_ASSERT(seq.size() > 0);
 
-	os << class_name() << "(";
-	seq[0].printraw(os);
-	os << ",indices=";
-	printrawindices(os);
-	os << ",hash=" << hashvalue << ",flags=" << flags << ")";
-}
-
-void indexed::printtree(std::ostream & os, unsigned indent) const
-{
-	debugmsg("indexed printtree", LOGLEVEL_PRINT);
-	GINAC_ASSERT(seq.size() > 0);
+	if (is_of_type(c, print_tree)) {
 
-	os << std::string(indent, ' ') << class_name() << ", " << seq.size()-1 << " indices";
-	os << ",hash=" << hashvalue << ",flags=" << flags << std::endl;
-	printtreeindices(os, indent);
-}
+		c.s << std::string(level, ' ') << class_name()
+		    << std::hex << ", hash=0x" << hashvalue << ", flags=0x" << flags << std::dec
+		    << ", " << seq.size()-1 << " indices";
+		switch (symmetry) {
+			case symmetric: c.s << ", symmetric"; break;
+			case antisymmetric: c.s << ", antisymmetric"; break;
+			default: break;
+		}
+		c.s << std::endl;
+		unsigned delta_indent = static_cast<const print_tree &>(c).delta_indent;
+		seq[0].print(c, level + delta_indent);
+		printindices(c, level + delta_indent);
 
-void indexed::print(std::ostream & os, unsigned upper_precedence) const
-{
-	debugmsg("indexed print", LOGLEVEL_PRINT);
-	GINAC_ASSERT(seq.size() > 0);
+	} else {
 
-	const ex & base = seq[0];
-	bool need_parens = is_ex_exactly_of_type(base, add) || is_ex_exactly_of_type(base, mul)
-	                || is_ex_exactly_of_type(base, ncmul) || is_ex_exactly_of_type(base, power);
-	if (need_parens)
-		os << "(";
-	os << base;
-	if (need_parens)
-		os << ")";
-	printindices(os);
+		const ex & base = seq[0];
+		bool need_parens = is_ex_exactly_of_type(base, add) || is_ex_exactly_of_type(base, mul)
+		                || is_ex_exactly_of_type(base, ncmul) || is_ex_exactly_of_type(base, power);
+		if (need_parens)
+			c.s << "(";
+		base.print(c);
+		if (need_parens)
+			c.s << ")";
+		printindices(c, level);
+	}
 }
 
 bool indexed::info(unsigned inf) const
@@ -360,6 +348,24 @@ ex indexed::eval(int level) const
 	return base.bp->eval_indexed(*this);
 }
 
+int indexed::degree(const ex & s) const
+{
+	return is_equal(*s.bp) ? 1 : 0;
+}
+
+int indexed::ldegree(const ex & s) const
+{
+	return is_equal(*s.bp) ? 1 : 0;
+}
+
+ex indexed::coeff(const ex & s, int n) const
+{
+	if (is_equal(*s.bp))
+		return n==1 ? _ex1() : _ex0();
+	else
+		return n==0 ? ex(*this) : _ex0();
+}
+
 ex indexed::thisexprseq(const exvector & v) const
 {
 	return indexed(symmetry, v);
@@ -400,38 +406,12 @@ ex indexed::expand(unsigned options) const
 // non-virtual functions in this class
 //////////
 
-void indexed::printrawindices(std::ostream & os) const
-{
-	if (seq.size() > 1) {
-		exvector::const_iterator it=seq.begin() + 1, itend = seq.end();
-		while (it != itend) {
-			it->printraw(os);
-			it++;
-			if (it != itend)
-				os << ",";
-		}
-	}
-}
-
-void indexed::printtreeindices(std::ostream & os, unsigned indent) const
-{
-	if (seq.size() > 1) {
-		exvector::const_iterator it=seq.begin() + 1, itend = seq.end();
-		while (it != itend) {
-			os << std::string(indent + delta_indent, ' ');
-			it->printraw(os);
-			os << std::endl;
-			it++;
-		}
-	}
-}
-
-void indexed::printindices(std::ostream & os) const
+void indexed::printindices(const print_context & c, unsigned level) const
 {
 	if (seq.size() > 1) {
 		exvector::const_iterator it=seq.begin() + 1, itend = seq.end();
 		while (it != itend) {
-			it->print(os);
+			it->print(c, level);
 			it++;
 		}
 	}
@@ -455,48 +435,6 @@ void indexed::assert_all_indices_of_type_idx(void) const
 // global functions
 //////////
 
-/** Given a vector of indices, split them into two vectors, one containing
- *  the free indices, the other containing the dummy indices. */
-static 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);
-}
-
 /** Check whether two sorted index vectors are consistent (i.e. equal). */
 static bool indices_consistent(const exvector & v1, const exvector & v2)
 {
@@ -515,6 +453,12 @@ static bool indices_consistent(const exvector & v1, const exvector & v2)
 	return true;
 }
 
+exvector indexed::get_indices(void) const
+{
+	GINAC_ASSERT(seq.size() >= 1);
+	return exvector(seq.begin() + 1, seq.end());
+}
+
 exvector indexed::get_dummy_indices(void) const
 {
 	exvector free_indices, dummy_indices;
@@ -522,6 +466,16 @@ exvector indexed::get_dummy_indices(void) const
 	return dummy_indices;
 }
 
+exvector indexed::get_dummy_indices(const indexed & other) const
+{
+	exvector indices = get_free_indices();
+	exvector other_indices = other.get_free_indices();
+	indices.insert(indices.end(), other_indices.begin(), other_indices.end());
+	exvector dummy_indices;
+	find_dummy_indices(indices, dummy_indices);
+	return dummy_indices;
+}
+
 exvector indexed::get_free_indices(void) const
 {
 	exvector free_indices, dummy_indices;
@@ -555,7 +509,7 @@ exvector mul::get_free_indices(void) const
 
 	// And remove the dummy indices
 	exvector free_indices, dummy_indices;
-	find_free_and_dummy(un.begin(), un.end(), free_indices, dummy_indices);
+	find_free_and_dummy(un, free_indices, dummy_indices);
 	return free_indices;
 }
 
@@ -570,7 +524,7 @@ exvector ncmul::get_free_indices(void) const
 
 	// And remove the dummy indices
 	exvector free_indices, dummy_indices;
-	find_free_and_dummy(un.begin(), un.end(), free_indices, dummy_indices);
+	find_free_and_dummy(un, free_indices, dummy_indices);
 	return free_indices;
 }
 
@@ -606,7 +560,7 @@ ex simplify_indexed_product(const ex & e, exvector & free_indices, const scalar_
 			} else if (is_ex_exactly_of_type(f, ncmul)) {
 				// Noncommutative factor found, split it as well
 				non_commutative = true; // everything becomes noncommutative, ncmul will sort out the commutative factors later
-				for (int j=0; j<f.nops(); i++)
+				for (int j=0; j<f.nops(); j++)
 					v.push_back(f.op(j));
 			} else
 				v.push_back(f);
@@ -623,33 +577,51 @@ try_again:
 		if (!is_ex_of_type(*it1, indexed))
 			continue;
 
-		// Indexed factor found, look for contraction candidates
+		// Indexed factor found, get free indices and look for contraction
+		// candidates
+		exvector free1, dummy1;
+		find_free_and_dummy(ex_to_indexed(*it1).seq.begin() + 1, ex_to_indexed(*it1).seq.end(), free1, dummy1);
+
 		exvector::iterator it2;
 		for (it2 = it1 + 1; it2 != itend; it2++) {
 
 			if (!is_ex_of_type(*it2, indexed))
 				continue;
 
+			// Find free indices of second factor and merge them with free
+			// indices of first factor
+			exvector un;
+			find_free_and_dummy(ex_to_indexed(*it2).seq.begin() + 1, ex_to_indexed(*it2).seq.end(), un, dummy1);
+			un.insert(un.end(), free1.begin(), free1.end());
+
 			// Check whether the two factors share dummy indices
-			exvector un(ex_to_indexed(*it1).seq.begin() + 1, ex_to_indexed(*it1).seq.end());
-			un.insert(un.end(), ex_to_indexed(*it2).seq.begin() + 1, ex_to_indexed(*it2).seq.end());
 			exvector free, dummy;
-			find_free_and_dummy(un.begin(), un.end(), free, dummy);
+			find_free_and_dummy(un, free, dummy);
 			if (dummy.size() == 0)
 				continue;
 
 			// At least one dummy index, is it a defined scalar product?
+			bool contracted = false;
 			if (free.size() == 0) {
 				if (sp.is_defined(*it1, *it2)) {
 					*it1 = sp.evaluate(*it1, *it2);
 					*it2 = _ex1();
-					something_changed = true;
-					goto try_again;
+					goto contraction_done;
 				}
 			}
 
+			// 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);
+			contracted = it1->op(0).bp->contract_with(it1, it2, v);
 			if (!contracted) {
 
 				// That didn't work; maybe the second object knows how to
@@ -657,11 +629,20 @@ try_again:
 				contracted = it2->op(0).bp->contract_with(it2, it1, v);
 			}
 			if (contracted) {
-				something_changed = true;
+contraction_done:
+				if (is_ex_exactly_of_type(*it1, add) || is_ex_exactly_of_type(*it2, add)
+				 || is_ex_exactly_of_type(*it1, mul) || is_ex_exactly_of_type(*it2, mul)) {
+
+					// One of the factors became a sum or product:
+					// re-expand expression and run again
+					ex r = non_commutative ? ex(ncmul(v)) : ex(mul(v));
+					return simplify_indexed(r, free_indices, sp);
+				}
 
 				// Both objects may have new indices now or they might
 				// even not be indexed objects any more, so we have to
 				// start over
+				something_changed = true;
 				goto try_again;
 			}
 		}
@@ -671,21 +652,24 @@ try_again:
 	exvector un, dummy_indices;
 	it1 = v.begin(); itend = v.end();
 	while (it1 != itend) {
-		if (is_ex_of_type(*it1, indexed)) {
-			const indexed & o = ex_to_indexed(*it1);
-			un.insert(un.end(), o.seq.begin() + 1, o.seq.end());
-		}
+		exvector free_indices_of_factor = it1->get_free_indices();
+		un.insert(un.end(), free_indices_of_factor.begin(), free_indices_of_factor.end());
 		it1++;
 	}
-	find_free_and_dummy(un.begin(), un.end(), free_indices, dummy_indices);
+	find_free_and_dummy(un, free_indices, dummy_indices);
 
-	if (something_changed) {
-		if (non_commutative)
-			return ncmul(v);
-		else
-			return mul(v);
-	} else
-		return e;
+	ex r;
+	if (something_changed)
+		r = non_commutative ? ex(ncmul(v)) : ex(mul(v));
+	else
+		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. */
@@ -705,15 +689,27 @@ ex simplify_indexed(const ex & e, exvector & free_indices, const scalar_products
 	// 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;