X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=blobdiff_plain;f=ginac%2Findexed.cpp;h=29dd86ae5e001f7c8903ce9861c702d95f90d9ec;hp=3fad73f8eb22ae2e533f8b889ac247dafbba5a2e;hb=172d1d7b2a12e2b3b6ff4ff35c37aeccafed3edb;hpb=f5e84af31b20c7f732bee375bacc152e7fb01e56

diff --git a/ginac/indexed.cpp b/ginac/indexed.cpp
index 3fad73f8..29dd86ae 100644
--- a/ginac/indexed.cpp
+++ b/ginac/indexed.cpp
@@ -3,7 +3,7 @@
  *  Implementation of GiNaC's indexed expressions. */
 
 /*
- *  GiNaC Copyright (C) 1999-2003 Johannes Gutenberg University Mainz, Germany
+ *  GiNaC Copyright (C) 1999-2004 Johannes Gutenberg University Mainz, Germany
  *
  *  This program is free software; you can redistribute it and/or modify
  *  it under the terms of the GNU General Public License as published by
@@ -36,6 +36,7 @@
 #include "lst.h"
 #include "archive.h"
 #include "utils.h"
+#include "integral.h"
 
 namespace GiNaC {
 
@@ -48,7 +49,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 +58,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 +106,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;
@@ -129,7 +130,7 @@ indexed::indexed(const symmetry & symm, const exvector & v, bool discardable) :
 	tinfo_key = TINFO_indexed;
 }
 
-indexed::indexed(const symmetry & symm, exvector * vp) : inherited(vp), symtree(symm)
+indexed::indexed(const symmetry & symm, std::auto_ptr<exvector> vp) : inherited(vp), symtree(symm)
 {
 	tinfo_key = TINFO_indexed;
 }
@@ -152,7 +153,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);
@@ -237,7 +238,7 @@ void indexed::do_print_latex(const print_latex & c, unsigned level) const
 
 void indexed::do_print_tree(const print_tree & c, unsigned level) const
 {
-	c.s << std::string(level, ' ') << class_name()
+	c.s << std::string(level, ' ') << class_name() << " @" << this
 	    << std::hex << ", hash=0x" << hashvalue << ", flags=0x" << flags << std::dec
 	    << ", " << seq.size()-1 << " indices"
 	    << ", symmetry=" << symtree << std::endl;
@@ -316,7 +317,7 @@ ex indexed::thiscontainer(const exvector & v) const
 	return indexed(ex_to<symmetry>(symtree), v);
 }
 
-ex indexed::thiscontainer(exvector * vp) const
+ex indexed::thiscontainer(std::auto_ptr<exvector> vp) const
 {
 	return indexed(ex_to<symmetry>(symtree), vp);
 }
@@ -325,20 +326,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 +500,34 @@ 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;
+}
+
+exvector integral::get_free_indices() const
+{
+	if (a.get_free_indices().size() || b.get_free_indices().size())
+		throw (std::runtime_error("integral::get_free_indices: boundary values should not have free indices"));
+	return f.get_free_indices();
 }
 
 /** Rename dummy indices in an expression.
@@ -1107,8 +1136,9 @@ 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() const
+ex ex::simplify_indexed(unsigned options) const
 {
 	exvector free_indices, dummy_indices;
 	scalar_products sp;
@@ -1121,8 +1151,9 @@ ex ex::simplify_indexed() 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) const
+ex ex::simplify_indexed(const scalar_products & sp, unsigned options) const
 {
 	exvector free_indices, dummy_indices;
 	return GiNaC::simplify_indexed(*this, free_indices, dummy_indices, sp);