Added integral class.

[ginac.git] / ginac / indexed.cpp

diff --git a/ginac/indexed.cpp b/ginac/indexed.cpp
index 562ca1d5c6c0bd95dbc0584df003b11cf1b505be..29dd86ae5e001f7c8903ce9861c702d95f90d9ec 100644 (file)
--- 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);
@@ -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);