]> www.ginac.de Git - ginac.git/blobdiff - ginac/indexed.cpp
* mul::expand() (mul.cpp): When multiplying two sums, be careful not to
[ginac.git] / ginac / indexed.cpp
index 5759b24d4c93e6f21e64a07e9ed392145e62db70..765a07c8bf4c7a1f8be7a17411a55dbcb64311b9 100644 (file)
@@ -3,7 +3,7 @@
  *  Implementation of GiNaC's indexed expressions. */
 
 /*
- *  GiNaC Copyright (C) 1999-2001 Johannes Gutenberg University Mainz, Germany
+ *  GiNaC Copyright (C) 1999-2002 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
@@ -20,8 +20,8 @@
  *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
  */
 
+#include <iostream>
 #include <stdexcept>
-#include <algorithm>
 
 #include "indexed.h"
 #include "idx.h"
 #include "print.h"
 #include "archive.h"
 #include "utils.h"
-#include "debugmsg.h"
 
 namespace GiNaC {
 
 GINAC_IMPLEMENT_REGISTERED_CLASS(indexed, exprseq)
 
 //////////
-// default constructor, destructor, copy constructor assignment operator and helpers
+// default ctor, dtor, copy ctor, assignment operator and helpers
 //////////
 
 indexed::indexed() : symtree(sy_none())
 {
-       debugmsg("indexed default constructor", LOGLEVEL_CONSTRUCT);
        tinfo_key = TINFO_indexed;
 }
 
@@ -64,63 +62,54 @@ DEFAULT_DESTROY(indexed)
 
 indexed::indexed(const ex & b) : inherited(b), symtree(sy_none())
 {
-       debugmsg("indexed constructor from ex", LOGLEVEL_CONSTRUCT);
        tinfo_key = TINFO_indexed;
        validate();
 }
 
 indexed::indexed(const ex & b, const ex & i1) : inherited(b, i1), symtree(sy_none())
 {
-       debugmsg("indexed constructor from ex,ex", LOGLEVEL_CONSTRUCT);
        tinfo_key = TINFO_indexed;
        validate();
 }
 
 indexed::indexed(const ex & b, const ex & i1, const ex & i2) : inherited(b, i1, i2), symtree(sy_none())
 {
-       debugmsg("indexed constructor from ex,ex,ex", LOGLEVEL_CONSTRUCT);
        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())
 {
-       debugmsg("indexed constructor from ex,ex,ex,ex", LOGLEVEL_CONSTRUCT);
        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())
 {
-       debugmsg("indexed constructor from ex,ex,ex,ex,ex", LOGLEVEL_CONSTRUCT);
        tinfo_key = TINFO_indexed;
        validate();
 }
 
 indexed::indexed(const ex & b, const symmetry & symm, const ex & i1, const ex & i2) : inherited(b, i1, i2), symtree(symm)
 {
-       debugmsg("indexed constructor from ex,symmetry,ex,ex", LOGLEVEL_CONSTRUCT);
        tinfo_key = TINFO_indexed;
        validate();
 }
 
 indexed::indexed(const ex & b, const symmetry & symm, const ex & i1, const ex & i2, const ex & i3) : inherited(b, i1, i2, i3), symtree(symm)
 {
-       debugmsg("indexed constructor from ex,symmetry,ex,ex,ex", LOGLEVEL_CONSTRUCT);
        tinfo_key = TINFO_indexed;
        validate();
 }
 
 indexed::indexed(const ex & b, const symmetry & symm, const ex & i1, const ex & i2, const ex & i3, const ex & i4) : inherited(b, i1, i2, i3, i4), symtree(symm)
 {
-       debugmsg("indexed constructor from ex,symmetry,ex,ex,ex,ex", LOGLEVEL_CONSTRUCT);
        tinfo_key = TINFO_indexed;
        validate();
 }
 
 indexed::indexed(const ex & b, const exvector & v) : inherited(b), symtree(sy_none())
 {
-       debugmsg("indexed constructor from ex,exvector", LOGLEVEL_CONSTRUCT);
        seq.insert(seq.end(), v.begin(), v.end());
        tinfo_key = TINFO_indexed;
        validate();
@@ -128,7 +117,6 @@ indexed::indexed(const ex & b, const exvector & v) : inherited(b), symtree(sy_no
 
 indexed::indexed(const ex & b, const symmetry & symm, const exvector & v) : inherited(b), symtree(symm)
 {
-       debugmsg("indexed constructor from ex,symmetry,exvector", LOGLEVEL_CONSTRUCT);
        seq.insert(seq.end(), v.begin(), v.end());
        tinfo_key = TINFO_indexed;
        validate();
@@ -136,19 +124,16 @@ indexed::indexed(const ex & b, const symmetry & symm, const exvector & v) : inhe
 
 indexed::indexed(const symmetry & symm, const exprseq & es) : inherited(es), symtree(symm)
 {
-       debugmsg("indexed constructor from symmetry,exprseq", LOGLEVEL_CONSTRUCT);
        tinfo_key = TINFO_indexed;
 }
 
 indexed::indexed(const symmetry & symm, const exvector & v, bool discardable) : inherited(v, discardable), symtree(symm)
 {
-       debugmsg("indexed constructor from symmetry,exvector", LOGLEVEL_CONSTRUCT);
        tinfo_key = TINFO_indexed;
 }
 
 indexed::indexed(const symmetry & symm, exvector * vp) : inherited(vp), symtree(symm)
 {
-       debugmsg("indexed constructor from symmetry,exvector *", LOGLEVEL_CONSTRUCT);
        tinfo_key = TINFO_indexed;
 }
 
@@ -158,7 +143,6 @@ indexed::indexed(const symmetry & symm, exvector * vp) : inherited(vp), symtree(
 
 indexed::indexed(const archive_node &n, const lst &sym_lst) : inherited(n, sym_lst)
 {
-       debugmsg("indexed constructor from archive_node", LOGLEVEL_CONSTRUCT);
        if (!n.find_ex("symmetry", symtree, sym_lst)) {
                // GiNaC versions <= 0.9.0 had an unsigned "symmetry" property
                unsigned symm = 0;
@@ -192,7 +176,6 @@ DEFAULT_UNARCHIVE(indexed)
 
 void indexed::print(const print_context & c, unsigned level) const
 {
-       debugmsg("indexed print", LOGLEVEL_PRINT);
        GINAC_ASSERT(seq.size() > 0);
 
        if (is_of_type(c, print_tree)) {
@@ -264,7 +247,7 @@ ex indexed::eval(int level) const
 
        // If the base object is 0, the whole object is 0
        if (base.is_zero())
-               return _ex0();
+               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)) {
@@ -282,7 +265,7 @@ ex indexed::eval(int level) const
                if (sig != INT_MAX) {
                        // Something has changed while sorting indices, more evaluations later
                        if (sig == 0)
-                               return _ex0();
+                               return _ex0;
                        return ex(sig) * thisexprseq(v);
                }
        }
@@ -291,24 +274,6 @@ ex indexed::eval(int level) const
        return ex_to<basic>(base).eval_indexed(*this);
 }
 
-int indexed::degree(const ex & s) const
-{
-       return is_equal(ex_to<basic>(s)) ? 1 : 0;
-}
-
-int indexed::ldegree(const ex & s) const
-{
-       return is_equal(ex_to<basic>(s)) ? 1 : 0;
-}
-
-ex indexed::coeff(const ex & s, int n) const
-{
-       if (is_equal(ex_to<basic>(s)))
-               return n==1 ? _ex1() : _ex0();
-       else
-               return n==0 ? ex(*this) : _ex0();
-}
-
 ex indexed::thisexprseq(const exvector & v) const
 {
        return indexed(ex_to<symmetry>(symtree), v);
@@ -327,7 +292,7 @@ ex indexed::expand(unsigned options) const
 
                // expand_indexed expands (a+b).i -> a.i + b.i
                const ex & base = seq[0];
-               ex sum = _ex0();
+               ex sum = _ex0;
                for (unsigned i=0; i<base.nops(); i++) {
                        exvector s = seq;
                        s[0] = base.op(i);
@@ -363,9 +328,10 @@ void indexed::printindices(const print_context & c, unsigned level) const
 
                        while (it != itend) {
                                bool cur_covariant = (is_ex_of_type(*it, varidx) ? ex_to<varidx>(*it).is_covariant() : true);
-                               if (first || cur_covariant != covariant) {
+                               if (first || cur_covariant != covariant) { // Variance changed
+                                       // The empty {} prevents indices from ending up on top of each other
                                        if (!first)
-                                               c.s << "}";
+                                               c.s << "}{}";
                                        covariant = cur_covariant;
                                        if (covariant)
                                                c.s << "_{";
@@ -415,7 +381,7 @@ void indexed::validate(void) const
  *  @see ex::diff */
 ex indexed::derivative(const symbol & s) const
 {
-       return _ex0();
+       return _ex0;
 }
 
 //////////
@@ -556,7 +522,6 @@ static ex rename_dummy_indices(const ex & e, exvector & global_dummy_indices, ex
                        }
                        it++;
                }
-               shaker_sort(global_dummy_indices.begin(), global_dummy_indices.end(), ex_is_less(), ex_swap());
 
                // If this is the first set of local indices, do nothing
                if (old_global_size == 0)
@@ -571,6 +536,7 @@ static ex rename_dummy_indices(const ex & e, exvector & global_dummy_indices, ex
        shaker_sort(local_syms.begin(), local_syms.end(), ex_is_less(), ex_swap());
        for (unsigned i=0; i<global_size; i++)
                global_syms.push_back(global_dummy_indices[i].op(0));
+       shaker_sort(global_syms.begin(), global_syms.end(), ex_is_less(), ex_swap());
 
        // Remove common indices
        exlist local_uniq, global_uniq;
@@ -601,13 +567,13 @@ ex simplify_indexed_product(const ex & e, exvector & free_indices, exvector & du
 
        if (is_ex_exactly_of_type(e, power)) {
                // We only get called for simple squares, split a^2 -> a*a
-               GINAC_ASSERT(e.op(1).is_equal(_ex2()));
+               GINAC_ASSERT(e.op(1).is_equal(_ex2));
                v.push_back(e.op(0));
                v.push_back(e.op(0));
        } else {
                for (unsigned i=0; i<e.nops(); i++) {
                        ex f = e.op(i);
-                       if (is_ex_exactly_of_type(f, power) && f.op(1).is_equal(_ex2())) {
+                       if (is_ex_exactly_of_type(f, power) && f.op(1).is_equal(_ex2)) {
                                v.push_back(f.op(0));
                    v.push_back(f.op(0));
                        } else if (is_ex_exactly_of_type(f, ncmul)) {
@@ -663,7 +629,7 @@ try_again:
                        if (free.empty()) {
                                if (sp.is_defined(*it1, *it2)) {
                                        *it1 = sp.evaluate(*it1, *it2);
-                                       *it2 = _ex1();
+                                       *it2 = _ex1;
                                        goto contraction_done;
                                }
                        }
@@ -735,7 +701,7 @@ contraction_done:
                        dummy_syms.append(local_dummy_indices[i].op(0));
                if (r.symmetrize(dummy_syms).is_zero()) {
                        free_indices.clear();
-                       return _ex0();
+                       return _ex0;
                }
        }
 
@@ -769,7 +735,7 @@ ex simplify_indexed(const ex & e, exvector & free_indices, exvector & dummy_indi
        // free indices in each term
        if (is_ex_exactly_of_type(e_expanded, add)) {
                bool first = true;
-               ex sum = _ex0();
+               ex sum = _ex0;
                free_indices.clear();
 
                for (unsigned i=0; i<e_expanded.nops(); i++) {
@@ -797,7 +763,7 @@ ex simplify_indexed(const ex & e, exvector & free_indices, exvector & dummy_indi
        // 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())))
+        || (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