- revamped indexed objects
authorChristian Bauer <Christian.Bauer@uni-mainz.de>
Tue, 6 Mar 2001 20:07:16 +0000 (20:07 +0000)
committerChristian Bauer <Christian.Bauer@uni-mainz.de>
Tue, 6 Mar 2001 20:07:16 +0000 (20:07 +0000)
- subs() works on matrices

34 files changed:
ginac/Makefile.am
ginac/add.cpp
ginac/add.h
ginac/basic.cpp
ginac/basic.h
ginac/constant.cpp
ginac/container.pl
ginac/ex.cpp
ginac/ex.h
ginac/expairseq.cpp
ginac/flags.h
ginac/function.pl
ginac/ginac.h
ginac/idx.cpp
ginac/idx.h
ginac/indexed.cpp
ginac/indexed.h
ginac/matrix.cpp
ginac/matrix.h
ginac/mul.cpp
ginac/mul.h
ginac/ncmul.cpp
ginac/ncmul.h
ginac/numeric.cpp
ginac/power.cpp
ginac/power.h
ginac/pseries.cpp
ginac/relational.cpp
ginac/structure.cpp
ginac/structure.pl
ginac/symbol.cpp
ginac/tensor.cpp [new file with mode: 0644]
ginac/tensor.h [new file with mode: 0644]
ginac/tinfos.h

index 7a3cd72..fffe951 100644 (file)
@@ -5,20 +5,17 @@ libginac_la_SOURCES = add.cpp archive.cpp basic.cpp constant.cpp ex.cpp \
   expairseq.cpp exprseq.cpp fail.cpp function.cpp inifcns.cpp \
   inifcns_trans.cpp inifcns_zeta.cpp inifcns_gamma.cpp matrix.cpp mul.cpp \
   normal.cpp numeric.cpp operators.cpp power.cpp registrar.cpp relational.cpp \
-  symbol.cpp pseries.cpp utils.cpp ncmul.cpp clifford.cpp structure.cpp \
-  color.cpp indexed.cpp idx.cpp isospin.cpp exprseq_suppl.cpp lst.cpp \
-  lst_suppl.cpp simp_lor.cpp coloridx.cpp lorentzidx.cpp lortensor.cpp \
-  input_parser.yy input_lexer.ll input_lexer.h \
-  remember.h remember.cpp debugmsg.h utils.h 
+  symbol.cpp pseries.cpp utils.cpp ncmul.cpp structure.cpp exprseq_suppl.cpp \
+  lst.cpp lst_suppl.cpp input_parser.yy input_lexer.ll input_lexer.h \
+  remember.h remember.cpp debugmsg.h utils.h idx.cpp indexed.cpp tensor.cpp
 libginac_la_LDFLAGS = -version-info $(LT_CURRENT):$(LT_REVISION):$(LT_AGE) \
   -release $(LT_RELEASE)
 ginacincludedir = $(includedir)/ginac
-ginacinclude_HEADERS = ginac.h add.h archive.h basic.h clifford.h color.h \
-  coloridx.h constant.h ex.h expair.h expairseq.h exprseq.h fail.h flags.h \
-  function.h idx.h indexed.h inifcns.h isospin.h lorentzidx.h lst.h matrix.h \
-  mul.h ncmul.h normal.h numeric.h operators.h power.h registrar.h \
-  relational.h pseries.h simp_lor.h structure.h symbol.h lortensor.h \
-  tinfos.h assertion.h version.h
+ginacinclude_HEADERS = ginac.h add.h archive.h basic.h constant.h ex.h \
+  expair.h expairseq.h exprseq.h fail.h flags.h function.h inifcns.h \
+  lst.h matrix.h mul.h ncmul.h normal.h numeric.h operators.h power.h \
+  registrar.h relational.h pseries.h structure.h symbol.h tinfos.h assertion.h \
+  version.h idx.h indexed.h tensor.h
 LFLAGS = -Pginac_yy -olex.yy.c
 YFLAGS = -p ginac_yy -d
 EXTRA_DIST = container.pl function.pl structure.pl input_parser.h version.h.in
index 0c22465..6c6a163 100644 (file)
@@ -354,17 +354,6 @@ ex add::eval(int level) const
        return this->hold();
 }
 
-exvector add::get_indices(void) const
-{
-       // FIXME: all terms in the sum should have the same indices (compatible
-       // tensors) however this is not checked, since there is no function yet
-       // which compares indices (idxvector can be unsorted)
-       if (seq.size()==0) {
-               return exvector();
-       }
-       return (seq.begin())->rest.get_indices();
-}    
-
 ex add::simplify_ncmul(const exvector & v) const
 {
        if (seq.size()==0) {
index 8a93186..8738bdf 100644 (file)
@@ -59,7 +59,7 @@ public:
        numeric integer_content(void) const;
        ex smod(const numeric &xi) const;
        numeric max_coefficient(void) const;
-       exvector get_indices(void) const;
+       exvector get_free_indices(void) const;
        ex simplify_ncmul(const exvector & v) const;
 protected:
        ex derivative(const symbol & s) const;
index 6adc4c4..2414c92 100644 (file)
@@ -117,7 +117,7 @@ void basic::archive(archive_node &n) const
 void basic::print(std::ostream & os, unsigned upper_precedence) const
 {
        debugmsg("basic print",LOGLEVEL_PRINT);
-       os << "[basic object]";
+       os << "[" << class_name() << " object]";
 }
 
 /** Output to ostream in ugly raw format, so brave developers can have a look
@@ -125,7 +125,7 @@ void basic::print(std::ostream & os, unsigned upper_precedence) const
 void basic::printraw(std::ostream & os) const
 {
        debugmsg("basic printraw",LOGLEVEL_PRINT);
-       os << "[basic object]";
+       os << "[" << class_name() << " object]";
 }
 
 /** Output to ostream formatted in tree- (indented-) form, so developers can
@@ -287,6 +287,31 @@ ex basic::evalf(int level) const
        return *this;
 }
 
+/** Perform automatic symbolic evaluations on indexed expression that
+ *  contains this object as the base expression. */
+ex basic::eval_indexed(const basic & i) const
+ // this function can't take a "const ex & i" because that would result
+ // in an infinite eval() loop
+{
+       // There is nothing to do for basic objects
+       return i.hold();
+}
+
+/** Try to contract two indexed expressions. If a contraction exists, the
+ *  function overwrites one or both arguments and returns true. Otherwise it
+ *  returns false. It is guaranteed that both expressions are of class
+ *  indexed (or a subclass) and that at least one dummy index has been
+ *  found.
+ *
+ *  @param self The first indexed expression; it's base object is *this
+ *  @param other The second indexed expression
+ *  @return true if the contraction was successful, false otherwise */
+bool basic::contract_with(ex & self, ex & other) const
+{
+       // Do nothing
+       return false;
+}
+
 /** Substitute a set of symbols by arbitrary expressions. The ex returned
  *  will already be evaluated. */
 ex basic::subs(const lst & ls, const lst & lr) const
@@ -320,7 +345,8 @@ ex basic::diff(const symbol & s, unsigned nth) const
        return ndiff;
 }
 
-exvector basic::get_indices(void) const
+/** Return a vector containing the free indices of an expression. */
+exvector basic::get_free_indices(void) const
 {
        return exvector(); // return an empty exvector
 }
index e132574..842417c 100644 (file)
@@ -127,8 +127,10 @@ public: // only const functions please (may break reference counting)
        virtual numeric integer_content(void) const;
        virtual ex smod(const numeric &xi) const;
        virtual numeric max_coefficient(void) const;
-       virtual exvector get_indices(void) const;
+       virtual exvector get_free_indices(void) const;
        virtual ex simplify_ncmul(const exvector & v) const;
+       virtual ex eval_indexed(const basic & i) const;
+       virtual bool contract_with(ex & self, ex & other) const;
 protected: // non-const functions should be called from class ex only
        virtual ex derivative(const symbol & s) const;
        virtual int compare_same_type(const basic & other) const;
index 56f5a85..9002515 100644 (file)
@@ -138,7 +138,7 @@ void constant::print(std::ostream & os, unsigned upper_precedence) const
 void constant::printraw(std::ostream & os) const
 {
        debugmsg("constant printraw",LOGLEVEL_PRINT);
-       os << "constant(" << name << ")";
+       os << class_name() << "(" << name << ")";
 }
 
 void constant::printtree(std::ostream & os, unsigned indent) const
index 80fb9fa..571330f 100755 (executable)
@@ -410,7 +410,7 @@ void ${CONTAINER}::printraw(std::ostream & os) const
 {
        debugmsg("${CONTAINER} printraw",LOGLEVEL_PRINT);
 
-       os << "${CONTAINER}(";
+       os << class_name() << "(";
        for (${STLT}::const_iterator cit=seq.begin(); cit!=seq.end(); ++cit) {
                (*cit).bp->printraw(os);
                os << ",";
index 6acde9f..5548255 100644 (file)
@@ -37,6 +37,7 @@
 #include "numeric.h"
 #include "power.h"
 #include "relational.h"
+#include "indexed.h"
 #include "input_lexer.h"
 #include "debugmsg.h"
 #include "utils.h"
@@ -241,10 +242,33 @@ ex ex::subs(const ex & e) const
        return bp->subs(e);
 }
 
-exvector ex::get_indices(void) const
+/** Return a vector containing the free indices of the object. */
+exvector ex::get_free_indices(void) const
 {
        GINAC_ASSERT(bp!=0);
-       return bp->get_indices();
+       return bp->get_free_indices();
+}
+
+/** 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
+{
+       return GiNaC::simplify_indexed(*this);
+}
+
+/** 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
+{
+       return GiNaC::simplify_indexed(*this, sp);
 }
 
 ex ex::simplify_ncmul(const exvector & v) const
index 0e2adaf..f2feba7 100644 (file)
@@ -35,6 +35,7 @@ extern const ex & _ex0(void);     ///<  single ex(numeric(0))
 
 class symbol;
 class lst;
+class scalar_products;
 
 /** Lightweight wrapper for GiNaC's symbolic objects.  Basically all it does is
  *  to hold a pointer to the other objects, manage the reference counting and
@@ -107,7 +108,9 @@ public:
        ex series(const ex & r, int order, unsigned options = 0) const;
        ex subs(const lst & ls, const lst & lr) const;
        ex subs(const ex & e) const;
-       exvector get_indices(void) const;
+       exvector get_free_indices(void) const;
+       ex simplify_indexed(void) const;
+       ex simplify_indexed(const scalar_products & sp) const;
        ex simplify_ncmul(const exvector & v) const;
        ex operator[](const ex & index) const;
        ex operator[](int i) const;
index bad422c..cfcfa16 100644 (file)
@@ -207,8 +207,7 @@ void expairseq::print(std::ostream &os, unsigned upper_precedence) const
 void expairseq::printraw(std::ostream &os) const
 {
        debugmsg("expairseq printraw",LOGLEVEL_PRINT);
-       
-       os << "expairseq(";
+       os << class_name() << "(";
        for (epvector::const_iterator cit=seq.begin(); cit!=seq.end(); ++cit) {
                os << "(";
                (*cit).rest.printraw(os);
index 93506f1..987a4df 100644 (file)
@@ -28,7 +28,8 @@ namespace GiNaC {
 class expand_options {
 public:
        enum {
-               expand_trigonometric = 0x0001
+               expand_trigonometric = 0x0001,
+               expand_indexed = 0x0002
        };
 };
 
@@ -128,13 +129,7 @@ public:
                has_indices,  // object has at least one index
 
                // answered by class idx
-               idx,
-
-               // answered by class coloridx
-               coloridx,
-
-               // answered by class lorentzidx
-               lorentzidx
+               idx
        };
 };
 
index 0a15766..ae9c832 100755 (executable)
@@ -642,7 +642,7 @@ void function::printraw(std::ostream & os) const
 
        GINAC_ASSERT(serial<registered_functions().size());
 
-       os << "function(name=" << registered_functions()[serial].name;
+       os << class_name() << "(name=" << registered_functions()[serial].name;
        for (exvector::const_iterator it=seq.begin(); it!=seq.end(); ++it) {
                os << ",";
                (*it).bp->print(os);
@@ -666,7 +666,7 @@ void function::printtree(std::ostream & os, unsigned indent) const
 
        GINAC_ASSERT(serial<registered_functions().size());
 
-       os << std::string(indent,' ') << "function "
+       os << std::string(indent,' ') << class_name() << " "
           << registered_functions()[serial].name
           << ", hash=" << hashvalue 
           << " (0x" << std::hex << hashvalue << std::dec << ")"
index b1a861a..5e08ce3 100644 (file)
@@ -33,7 +33,6 @@
 
 #include "constant.h"
 #include "fail.h"
-#include "idx.h"
 #include "lst.h"
 #include "matrix.h"
 #include "numeric.h"
 #include "operators.h"
 
 #ifndef GINAC_BASE_ONLY
+#include "idx.h"
 #include "indexed.h"
-#include "clifford.h"
-#include "coloridx.h"
-#include "color.h"
-#include "isospin.h"
-#include "lorentzidx.h"
-#include "simp_lor.h"
-#include "lortensor.h"
+#include "tensor.h"
 #endif // ndef GINAC_BASE_ONLY
 
 #ifdef __MAKECINT__
index 3c37884..89373cc 100644 (file)
@@ -23,9 +23,8 @@
 #include <stdexcept>
 
 #include "idx.h"
-#include "ex.h"
+#include "symbol.h"
 #include "lst.h"
-#include "relational.h"
 #include "archive.h"
 #include "utils.h"
 #include "debugmsg.h"
 namespace GiNaC {
 
 GINAC_IMPLEMENT_REGISTERED_CLASS(idx, basic)
+GINAC_IMPLEMENT_REGISTERED_CLASS(varidx, idx)
 
 //////////
 // default constructor, destructor, copy constructor assignment operator and helpers
 //////////
 
-// public
-
-idx::idx() : inherited(TINFO_idx), symbolic(true), covariant(false)
+idx::idx() : inherited(TINFO_idx)
 {
-       debugmsg("idx default constructor",LOGLEVEL_CONSTRUCT);
-       serial=next_serial++;
-       name=autoname_prefix()+ToString(serial);
+       debugmsg("idx default constructor", LOGLEVEL_CONSTRUCT);
 }
 
-// protected
+varidx::varidx() : covariant(false)
+{
+       debugmsg("varidx default constructor", LOGLEVEL_CONSTRUCT);
+       tinfo_key = TINFO_varidx;
+}
 
 void idx::copy(const idx & other)
 {
        inherited::copy(other);
-       serial=other.serial;
-       symbolic=other.symbolic;
-       name=other.name;
-       value=other.value;
-       covariant=other.covariant;
+       value = other.value;
+       dim = other.dim;
 }
 
-void idx::destroy(bool call_parent)
+void varidx::copy(const varidx & other)
 {
-       if (call_parent) inherited::destroy(call_parent);
+       inherited::copy(other);
+       covariant = other.covariant;
 }
 
-//////////
-// other constructors
-//////////
-
-// public
-
-/** Construct symbolic index, using an automatically generated unique name.
- *
- *  @param cov Index is covariant (contravariant otherwise)
- *  @return newly constructed index */
-idx::idx(bool cov) : inherited(TINFO_idx), symbolic(true), covariant(cov)
+void idx::destroy(bool call_parent)
 {
-       debugmsg("idx constructor from bool",LOGLEVEL_CONSTRUCT);
-       serial = next_serial++;
-       name = autoname_prefix()+ToString(serial);
+       if (call_parent)
+               inherited::destroy(call_parent);
 }
 
-/** Construct symbolic index with specified name.
- *
- *  @param n Symbolic index name
- *  @param cov Index is covariant (contravariant otherwise)
- *  @return newly constructed index */
-idx::idx(const std::string & n, bool cov) : inherited(TINFO_idx),  
-       symbolic(true), name(n), covariant(cov)
+void varidx::destroy(bool call_parent)
 {
-       debugmsg("idx constructor from string,bool",LOGLEVEL_CONSTRUCT);
-       serial = next_serial++;
+       if (call_parent)
+               inherited::destroy(call_parent);
 }
 
-/** Construct symbolic index with specified name.
- *
- *  @param n Symbolic index name
- *  @param cov Index is covariant (contravariant otherwise)
- *  @return newly constructed index */
-idx::idx(const char * n, bool cov) : inherited(TINFO_idx), symbolic(true), name(n), covariant(cov)
+//////////
+// other constructors
+//////////
+
+idx::idx(const ex & v, const ex & d) : inherited(TINFO_idx), value(v), dim(d)
 {
-       debugmsg("idx constructor from char*,bool",LOGLEVEL_CONSTRUCT);
-       serial = next_serial++;
+       debugmsg("idx constructor from ex,ex", LOGLEVEL_CONSTRUCT);
+       if (is_dim_numeric())
+               if (!dim.info(info_flags::posint))
+                       throw(std::invalid_argument("dimension of space must be a positive integer"));
 }
 
-/** Construct numeric index with specified value.
- *
- *  @param v Numeric index value
- *  @param cov Index is covariant (contravariant otherwise)
- *  @return newly constructed index */
-idx::idx(unsigned v, bool cov) : inherited(TINFO_idx), symbolic(false), value(v), covariant(cov)
+varidx::varidx(const ex & v, const ex & d, bool cov) : inherited(v, d), covariant(cov)
 {
-       debugmsg("idx constructor from unsigned,bool",LOGLEVEL_CONSTRUCT);
-       serial = 0;
+       debugmsg("varidx constructor from ex,ex,bool", LOGLEVEL_CONSTRUCT);
+       tinfo_key = TINFO_varidx;
 }
 
 //////////
 // archiving
 //////////
 
-/** Construct object from archive_node. */
 idx::idx(const archive_node &n, const lst &sym_lst) : inherited(n, sym_lst)
 {
        debugmsg("idx constructor from archive_node", LOGLEVEL_CONSTRUCT);
-       n.find_bool("symbolic", symbolic);
+       n.find_ex("value", value, sym_lst);
+       n.find_ex("dim", dim, sym_lst);
+}
+
+varidx::varidx(const archive_node &n, const lst &sym_lst) : inherited(n, sym_lst)
+{
+       debugmsg("varidx constructor from archive_node", LOGLEVEL_CONSTRUCT);
        n.find_bool("covariant", covariant);
-       if (symbolic) {
-               serial = next_serial++;
-               if (!(n.find_string("name", name)))
-                       name = autoname_prefix() + ToString(serial);
-       } else {
-               serial = 0;
-               n.find_unsigned("value", value);
-       }
 }
 
-/** Unarchive the object. */
 ex idx::unarchive(const archive_node &n, const lst &sym_lst)
 {
-       ex s = (new idx(n, sym_lst))->setflag(status_flags::dynallocated);
+       return (new idx(n, sym_lst))->setflag(status_flags::dynallocated);
+}
 
-       if (ex_to_idx(s).symbolic) {
-               // If idx is in sym_lst, return the existing idx
-               for (unsigned i=0; i<sym_lst.nops(); i++) {
-                       if (is_ex_of_type(sym_lst.op(i), idx) && (ex_to_idx(sym_lst.op(i)).name == ex_to_idx(s).name))
-                               return sym_lst.op(i);
-               }
-       }
-       return s;
+ex varidx::unarchive(const archive_node &n, const lst &sym_lst)
+{
+       return (new varidx(n, sym_lst))->setflag(status_flags::dynallocated);
 }
 
-/** Archive the object. */
 void idx::archive(archive_node &n) const
 {
        inherited::archive(n);
-       n.add_bool("symbolic", symbolic);
+       n.add_ex("value", value);
+       n.add_ex("dim", dim);
+}
+
+void varidx::archive(archive_node &n) const
+{
+       inherited::archive(n);
        n.add_bool("covariant", covariant);
-       if (symbolic)
-               n.add_string("name", name);
-       else
-               n.add_unsigned("value", value);
 }
 
 //////////
 // functions overriding virtual functions from bases classes
 //////////
 
-// public
-
 void idx::printraw(std::ostream & os) const
 {
-       debugmsg("idx printraw",LOGLEVEL_PRINT);
-
-       os << "idx(";
-
-       if (symbolic) {
-               os << "symbolic,name=" << name;
-       } else {
-               os << "non symbolic,value=" << value;
-       }
+       debugmsg("idx printraw", LOGLEVEL_PRINT);
 
-       if (covariant) {
-               os << ",covariant";
-       } else {
-               os << ",contravariant";
-       }
-
-       os << ",serial=" << serial;
+       os << class_name() << "(";
+       value.printraw(os);
+       os << ",dim=";
+       dim.printraw(os);
        os << ",hash=" << hashvalue << ",flags=" << flags;
        os << ")";
 }
@@ -195,334 +152,173 @@ void idx::printtree(std::ostream & os, unsigned indent) const
 {
        debugmsg("idx printtree",LOGLEVEL_PRINT);
 
-       os << std::string(indent,' ') << "idx: ";
-
-       if (symbolic) {
-               os << "symbolic,name=" << name;
-       } else {
-               os << "non symbolic,value=" << value;
-       }
-
-       if (covariant) {
-               os << ",covariant";
-       } else {
-               os << ",contravariant";
-       }
-
-       os << ", serial=" << serial
-          << ", hash=" << hashvalue
+       os << std::string(indent, ' ') << "type=" << class_name();
+       value.printtree(os, indent + delta_indent);
+       os << std::string(indent, ' ');
+       os << ", hash=" << hashvalue
           << " (0x" << std::hex << hashvalue << std::dec << ")"
           << ", flags=" << flags << std::endl;
 }
 
 void idx::print(std::ostream & os, unsigned upper_precedence) const
 {
-       debugmsg("idx print",LOGLEVEL_PRINT);
+       debugmsg("idx print", LOGLEVEL_PRINT);
 
-       if (covariant) {
-               os << "_";
-       } else {
-               os << "~";
-       }
-       if (symbolic) {
-               os << name;
-       } else {
-               os << value;
-       }
-}
+       os << ".";
 
-bool idx::info(unsigned inf) const
-{
-       if (inf==info_flags::idx) return true;
-       return inherited::info(inf);
+       bool need_parens = !(is_ex_exactly_of_type(value, numeric) || is_ex_of_type(value, symbol));
+       if (need_parens)
+               os << "(";
+       os << value;
+       if (need_parens)
+               os << ")";
 }
 
-ex idx::subs(const lst & ls, const lst & lr) const
+void varidx::print(std::ostream & os, unsigned upper_precedence) const
 {
-       GINAC_ASSERT(ls.nops()==lr.nops());
-#ifdef DO_GINAC_ASSERT
-       for (unsigned i=0; i<ls.nops(); i++) {
-               GINAC_ASSERT(is_ex_exactly_of_type(ls.op(i),symbol)||
-                          is_ex_of_type(ls.op(i),idx));
-       }
-#endif // def DO_GINAC_ASSERT
+       debugmsg("varidx print", LOGLEVEL_PRINT);
 
-       for (unsigned i=0; i<ls.nops(); i++) {
-               if (is_equal(*(ls.op(i)).bp)) {
-                       return lr.op(i);
-               }
-       }
-       return *this;
+       if (covariant)
+               os << ".";
+       else
+               os << "~";
+
+       bool need_parens = !(is_ex_exactly_of_type(value, numeric) || is_ex_of_type(value, symbol));
+       if (need_parens)
+               os << "(";
+       os << value;
+       if (need_parens)
+               os << ")";
 }
 
-// protected
+bool idx::info(unsigned inf) const
+{
+       if (inf == info_flags::idx)
+               return true;
+       return inherited::info(inf);
+}
 
+/** Returns order relation between two indices of the same type. The order
+ *  must be such that dummy indices lie next to each other. */
 int idx::compare_same_type(const basic & other) const
 {
-       GINAC_ASSERT(is_of_type(other,idx));
+       GINAC_ASSERT(is_of_type(other, idx));
        const idx &o = static_cast<const idx &>(other);
 
-       if (covariant!=o.covariant) {
-               // different co/contravariant
-               return covariant ? -1 : 1;
-       }
+       int cmpval = value.compare(o.value);
+       if (cmpval)
+               return cmpval;
+       return dim.compare(o.dim);
+}
 
-       if ((!symbolic) && (!o.symbolic)) {
-               // non-symbolic, of equal type: compare values
-               if (value==o.value) {
-                       return 0;
-               }
-               return value<o.value ? -1 : 1;
-       }
+int varidx::compare_same_type(const basic & other) const
+{
+       GINAC_ASSERT(is_of_type(other, varidx));
+       const varidx &o = static_cast<const varidx &>(other);
 
-       if (symbolic && o.symbolic) {
-               // both symbolic: compare serials
-               if (serial==o.serial) {
-                       return 0;
-               }
-               return serial<o.serial ? -1 : 1;
-       }
+       int cmpval = inherited::compare_same_type(other);
+       if (cmpval)
+               return cmpval;
 
-       // one symbolic, one value: value is sorted first
-       return o.symbolic ? -1 : 1;
+       // Check variance last so dummy indices will end up next to each other
+       if (covariant != o.covariant)
+               return covariant ? -1 : 1;
+       return 0;
 }
 
-bool idx::is_equal_same_type(const basic & other) const
+ex idx::subs(const lst & ls, const lst & lr) const
 {
-       GINAC_ASSERT(is_of_type(other,idx));
-       const idx &o = static_cast<const idx &>(other);
+       GINAC_ASSERT(ls.nops() == lr.nops());
 
-       if (covariant != o.covariant) return false;
-       if (symbolic != o.symbolic) return false;
-       if (symbolic && o.symbolic) return serial==o.serial;
-       return value==o.value;
-}    
+       // First look for index substitutions
+       for (unsigned i=0; i<ls.nops(); i++) {
+               if (is_equal(*(ls.op(i)).bp)) {
 
-unsigned idx::calchash(void) const
-{
-       hashvalue=golden_ratio_hash(golden_ratio_hash(tinfo_key ^ serial));
-       setflag(status_flags::hash_calculated);
-       return hashvalue;
-}
+                       // Substitution index->index
+                       if (is_ex_of_type(lr.op(i), idx))
+                               return lr.op(i);
 
-//////////
-// new virtual functions which can be overridden by derived classes
-//////////
+                       // Otherwise substitute value
+                       idx *i_copy = static_cast<idx *>(duplicate());
+                       i_copy->value = lr.op(i);
+                       return i_copy->setflag(status_flags::dynallocated);
+               }
+       }
 
-// public
+       // None, substitute objects in value (not in dimension)
+       const ex &subsed_value = value.subs(ls, lr);
+       if (are_ex_trivially_equal(value, subsed_value))
+               return *this;
 
-/** Check whether the index forms a co-/contravariant pair with another
- *  index (i.e. same name/value but opposite co-/contravariance). */
-bool idx::is_co_contra_pair(const basic & other) const
-{
-       // like is_equal_same_type(), but tests for different covariant status
-       GINAC_ASSERT(is_of_type(other,idx));
-       const idx & o=static_cast<const idx &>(const_cast<basic &>(other));
-
-       if (covariant==o.covariant) return false;
-       if (symbolic!=o.symbolic) return false;
-       if (symbolic && o.symbolic) return serial==o.serial;
-       return value==o.value;
-}    
-
-/** Toggle co-/contravariance of index. */
-ex idx::toggle_covariant(void) const
-{
-       idx * i_copy=static_cast<idx *>(duplicate());
-       i_copy->covariant = !i_copy->covariant;
-       i_copy->clearflag(status_flags::hash_calculated);
+       idx *i_copy = static_cast<idx *>(duplicate());
+       i_copy->value = subsed_value;
        return i_copy->setflag(status_flags::dynallocated);
 }
 
 //////////
-// non-virtual functions in this class
+// new virtual functions
 //////////
 
-// private
-
-std::string & idx::autoname_prefix(void)
+bool idx::is_dummy_pair_same_type(const basic & other) const
 {
-       static std::string * s = new std::string("index");
-       return *s;
-}
-
-//////////
-// static member variables
-//////////
+       const idx &o = static_cast<const idx &>(other);
 
-// protected
+       // Only pure symbols form dummy pairs, "2n+1" doesn't
+       if (!is_ex_of_type(value, symbol))
+               return false;
 
-unsigned idx::next_serial=0;
+       // Value must be equal, of course
+       if (!value.is_equal(o.value))
+               return false;
 
-//////////
-// other functions
-//////////
+       // Also the dimension
+       return dim.is_equal(o.dim);
+}
 
-/** Bring a vector of indices into a canonic order. This operation only makes
- *  sense if the object carrying these indices is either symmetric or totally
- *  antisymmetric with respect to the indices.
- *
- *  @param iv Index vector
- *  @param antisymmetric Whether the object carrying the indices is antisymmetric (symmetric otherwise)
- *  @return the sign introduced by the reordering of the indices. For symmetric
- *          objects this is always +1. For antisymmetric objects this is either
- *          +1 or -1 or 0 (if two equal indices were encountered). If the index
- *          vector was unchanged this function returns INT_MAX. */
-int canonicalize_indices(exvector & iv, bool antisymmetric)
+bool varidx::is_dummy_pair_same_type(const basic & other) const
 {
-       if (iv.size()<2) {
-               // nothing do to for 0 or 1 indices
-               return INT_MAX;
-       }
+       const varidx &o = static_cast<const varidx &>(other);
 
-       bool something_changed=false;
-       int sig=1;
-
-       // simple bubble sort algorithm should be sufficient for the small number of indices needed
-       exvector::const_iterator last_idx=iv.end();
-       exvector::const_iterator next_to_last_idx=iv.end()-1;
-       for (exvector::iterator it1=iv.begin(); it1!=next_to_last_idx; ++it1) {
-               for (exvector::iterator it2=it1+1; it2!=last_idx; ++it2) {
-                       int cmpval=(*it1).compare(*it2);
-                       if (cmpval==1) {
-                               iter_swap(it1,it2);
-                               something_changed=true;
-                               if (antisymmetric) sig=-sig;
-                       } else if ((cmpval==0) && antisymmetric) {
-                               something_changed=true;
-                               sig=0;
-                       }
-               }
-       }
+       // Variance must be opposite
+       if (covariant == o.covariant)
+               return false;
 
-       return something_changed ? sig : INT_MAX;
+       return inherited::is_dummy_pair_same_type(other);
 }
 
-/** Build a vector of indices as the set intersection of two other index
- *  vectors (i.e. the returned vector contains the indices which appear in
- *  both source vectors). */
-exvector idx_intersect(const exvector & iv1, const exvector & iv2)
-{
-       // Create union vector
-       exvector iv_union;
-       iv_union.reserve(iv1.size() + iv2.size());
-       iv_union.insert(iv_union.end(), iv1.begin(), iv1.end());
-       iv_union.insert(iv_union.end(), iv2.begin(), iv2.end());
-
-       // Sort it
-       canonicalize_indices(iv_union);
-
-       // Look for duplicates
-       exvector iv_intersect;
-       exvector::const_iterator cit = iv_union.begin(), citend = iv_union.end();
-       ex e;
-       if (cit != citend)
-               e = *cit++;
-       while (cit != citend) {
-               if (e.is_equal(*cit)) {
-                       iv_intersect.push_back(e);
-                       do {
-                               cit++;
-                       } while (cit != citend && e.is_equal(*cit));
-                       if (cit == citend)
-                               break;
-               }
-               e = *cit++;
-       }
-       return iv_intersect;
-}
+//////////
+// non-virtual functions
+//////////
 
-/** Given a vector iv3 of three indices and a vector iv2 of two indices
- *  where iv2 is a subset of iv3, return the (free) index that is in iv3
- *  but not in iv2 and the sign introduced by permuting that index to the
- *  front.
- *
- *  @param iv3 Vector of 3 indices
- *  @param iv2 Vector of 2 indices, must be a subset of iv3
- *  @param sig Returns the sign introduced by permuting the free index to the
- *             front if the object carrying the indices was antisymmetric (if
- *             it's symmetric, you can just ignore the returned value).
- *  @return the free index (the one that is in iv3 but not in iv2) */
-ex permute_free_index_to_front(const exvector & iv3, const exvector & iv2, int * sig)
+ex varidx::toggle_variance(void) const
 {
-       // match (return value,iv2) to iv3 by permuting indices
-       // iv3 is always cyclic
-
-       GINAC_ASSERT(iv3.size()==3);
-       GINAC_ASSERT(iv2.size()==2);
+       varidx *i_copy = static_cast<varidx *>(duplicate());
+       i_copy->covariant = !i_copy->covariant;
+       i_copy->clearflag(status_flags::hash_calculated);
+       return i_copy->setflag(status_flags::dynallocated);
+}
 
-       *sig=1;
+//////////
+// global functions
+//////////
 
-#define TEST_PERMUTATION(A,B,C,P) \
-       if ((iv3[B].is_equal(iv2[0]))&&(iv3[C].is_equal(iv2[1]))) { \
-               *sig=P; \
-               return iv3[A]; \
-       }
-       
-       TEST_PERMUTATION(0,1,2,  1);
-       TEST_PERMUTATION(0,2,1, -1);
-       TEST_PERMUTATION(1,0,2, -1);
-       TEST_PERMUTATION(1,2,0,  1);
-       TEST_PERMUTATION(2,0,1,  1);
-       TEST_PERMUTATION(2,1,0, -1);
-       throw(std::logic_error("permute_free_index_to_front(): no valid permutation found"));
-}
-       
-/** Substitute one index in a vector of expressions.
- *
- *  @param v Vector to substitute in (will be modified)
- *  @param is Index being substituted
- *  @param ir Index to replace by
- *  @return number of performed substitutions */
-unsigned subs_index_in_exvector(exvector & v, const ex & is, const ex & ir)
+bool is_dummy_pair(const idx & i1, const idx & i2)
 {
-       exvector::iterator it;
-       unsigned replacements=0;
-       unsigned current_replacements;
-
-       GINAC_ASSERT(is_ex_of_type(is,idx));
-       GINAC_ASSERT(is_ex_of_type(ir,idx));
-   
-       for (it=v.begin(); it!=v.end(); ++it) {
-               current_replacements=count_index(*it,is);
-               if (current_replacements>0) {
-                       (*it)=(*it).subs(is==ir);
-               }
-               replacements += current_replacements;
-       }
-       return replacements;
-}
+       // The indices must be of exactly the same type
+       if (i1.tinfo() != i2.tinfo())
+               return false;
 
-/** Count number of times a given index appears in the index vector of an
- *  indexed object.
- *
- *  @param e Indexed object
- *  @param i Index to look for
- *  @return number of times the index was found */
-unsigned count_index(const ex & e, const ex & i)
-{
-       exvector idxv=e.get_indices();
-       unsigned count=0;
-       for (exvector::const_iterator cit=idxv.begin(); cit!=idxv.end(); ++cit) {
-               if ((*cit).is_equal(i)) count++;
-       }
-       return count;
+       // Same type, let the indices decide whether they are paired
+       return i1.is_dummy_pair_same_type(i2);
 }
 
-/** Substitute multiple indices in an expression.
- *
- *  @param e Expression to substitute in
- *  @param idxv_subs Vector of indices being substituted
- *  @param idxv_repl Vector of indices to replace by (1:1 correspondence to idxv_subs)
- *  @return expression with substituted indices */
-ex subs_indices(const ex & e, const exvector & idxv_subs, const exvector & idxv_repl)
+bool is_dummy_pair(const ex & e1, const ex & e2)
 {
-       GINAC_ASSERT(idxv_subs.size()==idxv_repl.size());
-       ex res=e;
-       for (unsigned i=0; i<idxv_subs.size(); ++i) {
-               res=res.subs(idxv_subs[i]==idxv_repl[i]);
-       }
-       return res;
+       // The expressions must be indices
+       if (!is_ex_of_type(e1, idx) || !is_ex_of_type(e2, idx))
+               return false;
+
+       return is_dummy_pair(ex_to_idx(e1), ex_to_idx(e2));
 }
 
 } // namespace GiNaC
index 66e2502..c760f67 100644 (file)
 #ifndef __GINAC_IDX_H__
 #define __GINAC_IDX_H__
 
-#include <string>
-//#include <vector>
-#include "basic.h"
 #include "ex.h"
 
 namespace GiNaC {
 
 
-/** This class holds one index of an indexed object. Indices can be symbolic
- *  (e.g. "mu", "i") or numeric (unsigned integer), and they can be contravariant
- *  (the default) or covariant. */
+/** This class holds one index of an indexed object. Indices can
+ *  theoretically consist of any symbolic expression but they are usually
+ *  only just a symbol (e.g. "mu", "i") or numeric (integer). Indices belong
+ *  to a space with a certain numeric or symbolic dimension. */
 class idx : public basic
 {
        GINAC_DECLARE_REGISTERED_CLASS(idx, basic)
 
        // other constructors
 public:
-       explicit idx(bool cov);
-       explicit idx(const std::string & n, bool cov=false);
-       explicit idx(const char * n, bool cov=false);
-       explicit idx(unsigned v, bool cov=false); 
+       /** Construct index with given value and dimension.
+        *
+        *  @param v Value of index (numeric or symbolic)
+        *  @param dim Dimension of index space (numeric or symbolic)
+        *  @return newly constructed index */
+       explicit idx(const ex & v, const ex & dim);
 
        // functions overriding virtual functions from bases classes
 public:
@@ -52,56 +52,96 @@ public:
        void print(std::ostream & os, unsigned upper_precedence=0) const;
        bool info(unsigned inf) const;
 protected:
-       bool is_equal_same_type(const basic & other) const;
-       unsigned calchash(void) const;
        ex subs(const lst & ls, const lst & lr) const;
 
-       // new virtual functions which can be overridden by derived classes
+       // new virtual functions in this class
 public:
-       virtual bool is_co_contra_pair(const basic & other) const;
-       virtual ex toggle_covariant(void) const;
+       /** Check whether the index forms a dummy index pair with another index
+        *  of the same type. */
+       virtual bool is_dummy_pair_same_type(const basic & other) const;
 
        // non-virtual functions in this class
 public:
-       /** Check whether index is symbolic (not numeric). */
-       bool is_symbolic(void) const {return symbolic;}
+       /** Get value of index. */
+       ex get_value(void) const {return value;}
 
-       /** Get numeric value of index. Undefined for symbolic indices. */
-       unsigned get_value(void) const {return value;}
+       /** Check whether the index is numeric. */
+       bool is_numeric(void) const {return is_ex_exactly_of_type(value, numeric);}
 
-       /** Check whether index is covariant (not contravariant). */
+       /** Check whether the index is symbolic. */
+       bool is_symbolic(void) const {return !is_ex_exactly_of_type(value, numeric);}
+
+       /** Get dimension of index space. */
+       ex get_dim(void) const {return dim;}
+
+       /** Check whether the dimension is numeric. */
+       bool is_dim_numeric(void) const {return is_ex_exactly_of_type(dim, numeric);}
+
+       /** Check whether the dimension is symbolic. */
+       bool is_dim_symbolic(void) const {return !is_ex_exactly_of_type(dim, numeric);}
+
+       // member variables
+protected:
+       ex value; /**< Expression that constitutes the index (numeric or symbolic name) */
+       ex dim;   /**< Dimension of space (can be symbolic or numeric) */
+};
+
+
+/** This class holds an index with a variance (co- or contravariant). There
+ *  is an associated metric tensor that can be used to raise/lower indices. */
+class varidx : public idx
+{
+       GINAC_DECLARE_REGISTERED_CLASS(varidx, idx)
+
+       // other constructors
+public:
+       /** Construct index with given value, dimension and variance.
+        *
+        *  @param v Value of index (numeric or symbolic)
+        *  @param dim Dimension of index space (numeric or symbolic)
+        *  @param covariant Make covariant index (default is contravariant)
+        *  @return newly constructed index */
+       varidx(const ex & v, const ex & dim, bool covariant = false);
+
+       // functions overriding virtual functions from bases classes
+public:
+       void print(std::ostream & os, unsigned upper_precedence=0) const;
+       bool is_dummy_pair_same_type(const basic & other) const;
+
+       // non-virtual functions in this class
+public:
+       /** Check whether the index is covariant. */
        bool is_covariant(void) const {return covariant;}
 
-       void setname(const std::string & n) {name=n;}
-       std::string getname(void) const {return name;}
+       /** Check whether the index is contravariant (not covariant). */
+       bool is_contravariant(void) const {return !covariant;}
 
-private:
-       std::string & autoname_prefix(void);
+       /** Make a new index with the same value but the opposite variance. */
+       ex toggle_variance(void) const;
 
        // member variables
 protected:
-       unsigned serial;  /**< Unique serial number for comparing symbolic indices */
-       bool symbolic;    /**< Is index symbolic? */
-       std::string name; /**< Symbolic name (if symbolic == true) */
-       unsigned value;   /**< Numeric value (if symbolic == false) */
-       static unsigned next_serial;
-       bool covariant;   /**< x_mu, default is contravariant: x~mu */
+       bool covariant; /**< x.mu, default is contravariant: x~mu */
 };
 
+
 // utility functions
-inline const idx &ex_to_idx(const ex &e)
+inline const idx &ex_to_idx(const ex & e)
 {
        return static_cast<const idx &>(*e.bp);
 }
 
-// global functions
+inline const varidx &ex_to_varidx(const ex & e)
+{
+       return static_cast<const varidx &>(*e.bp);
+}
+
+/** Check whether two indices form a dummy pair. */
+bool is_dummy_pair(const idx & i1, const idx & i2);
+
+/** Check whether two expressions form a dummy index pair. */
+bool is_dummy_pair(const ex & e1, const ex & e2);
 
-int canonicalize_indices(exvector & iv, bool antisymmetric=false);
-exvector idx_intersect(const exvector & iv1, const exvector & iv2);
-ex permute_free_index_to_front(const exvector & iv3, const exvector & iv2, int * sig);
-unsigned subs_index_in_exvector(exvector & v, const ex & is, const ex & ir);
-ex subs_indices(const ex & e, const exvector & idxv_contra, const exvector & idxv_co);
-unsigned count_index(const ex & e, const ex & i);
 
 } // namespace GiNaC
 
index 9456f38..ea4266c 100644 (file)
@@ -1,6 +1,6 @@
 /** @file indexed.cpp
  *
- *  Implementation of GiNaC's index carrying objects. */
+ *  Implementation of GiNaC's indexed expressions. */
 
 /*
  *  GiNaC Copyright (C) 1999-2001 Johannes Gutenberg University Mainz, Germany
  *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
  */
 
-#include <string>
+#include <stdexcept>
 
 #include "indexed.h"
-#include "ex.h"
 #include "idx.h"
+#include "add.h"
+#include "mul.h"
+#include "ncmul.h"
+#include "power.h"
+#include "archive.h"
+#include "utils.h"
 #include "debugmsg.h"
 
 namespace GiNaC {
@@ -35,107 +40,119 @@ GINAC_IMPLEMENT_REGISTERED_CLASS(indexed, exprseq)
 // default constructor, destructor, copy constructor assignment operator and helpers
 //////////
 
-// public
-
-indexed::indexed()
+indexed::indexed() : symmetry(unknown)
 {
-       debugmsg("indexed default constructor",LOGLEVEL_CONSTRUCT);
-       tinfo_key=TINFO_indexed;
+       debugmsg("indexed default constructor", LOGLEVEL_CONSTRUCT);
+       tinfo_key = TINFO_indexed;
 }
 
-// protected
-
 void indexed::copy(const indexed & other)
 {
        inherited::copy(other);
+       symmetry = other.symmetry;
 }
 
 void indexed::destroy(bool call_parent)
 {
-       if (call_parent) {
+       if (call_parent)
                inherited::destroy(call_parent);
-       }
 }
 
 //////////
 // other constructors
 //////////
 
-// public
+indexed::indexed(const ex & b) : inherited(b), symmetry(unknown)
+{
+       debugmsg("indexed constructor from ex", LOGLEVEL_CONSTRUCT);
+       tinfo_key = TINFO_indexed;
+       GINAC_ASSERT(all_indices_of_type_idx());
+}
 
-/** Construct indexed object with one index. The index must be of class idx
- *  or a subclass.
- *
- *  @param i1 The index
- *  @return newly constructed indexed object */
-indexed::indexed(const ex & i1) : inherited(i1)
+indexed::indexed(const ex & b, const ex & i1) : inherited(b, i1), symmetry(unknown)
 {
-       debugmsg("indexed constructor from ex",LOGLEVEL_CONSTRUCT);
-       tinfo_key=TINFO_indexed;
-       GINAC_ASSERT(all_of_type_idx());
+       debugmsg("indexed constructor from ex,ex", LOGLEVEL_CONSTRUCT);
+       tinfo_key = TINFO_indexed;
+       GINAC_ASSERT(all_indices_of_type_idx());
 }
 
-/** Construct indexed object with two indices. The indices must be of class
- *  idx or a subclass.
- *
- *  @param i1 First index
- *  @param i2 Second index
- *  @return newly constructed indexed object */
-indexed::indexed(const ex & i1, const ex & i2) : inherited(i1,i2)
+indexed::indexed(const ex & b, const ex & i1, const ex & i2) : inherited(b, i1, i2), symmetry(unknown)
 {
-       debugmsg("indexed constructor from ex,ex",LOGLEVEL_CONSTRUCT);
-       tinfo_key=TINFO_indexed;
-       GINAC_ASSERT(all_of_type_idx());
+       debugmsg("indexed constructor from ex,ex,ex", LOGLEVEL_CONSTRUCT);
+       tinfo_key = TINFO_indexed;
+       GINAC_ASSERT(all_indices_of_type_idx());
 }
 
-/** Construct indexed object with three indices. The indices must be of class
- *  idx or a subclass.
- *
- *  @param i1 First index
- *  @param i2 Second index
- *  @param i3 Third index
- *  @return newly constructed indexed object */
-indexed::indexed(const ex & i1, const ex & i2, const ex & i3)
-  : inherited(i1,i2,i3)
+indexed::indexed(const ex & b, const ex & i1, const ex & i2, const ex & i3) : inherited(b, i1, i2, i3), symmetry(unknown)
 {
-       debugmsg("indexed constructor from ex,ex,ex",LOGLEVEL_CONSTRUCT);
-       tinfo_key=TINFO_indexed;
-       GINAC_ASSERT(all_of_type_idx());
+       debugmsg("indexed constructor from ex,ex,ex,ex", LOGLEVEL_CONSTRUCT);
+       tinfo_key = TINFO_indexed;
+       GINAC_ASSERT(all_indices_of_type_idx());
 }
 
-/** Construct indexed object with four indices. The indices must be of class
- *  idx or a subclass.
- *
- *  @param i1 First index
- *  @param i2 Second index
- *  @param i3 Third index
- *  @param i4 Fourth index
- *  @return newly constructed indexed object */
-indexed::indexed(const ex & i1, const ex & i2, const ex & i3, const ex & i4)
-  : inherited(i1,i2,i3,i4)
+indexed::indexed(const ex & b, const ex & i1, const ex & i2, const ex & i3, const ex & i4) : inherited(b, i1, i2, i3, i4), symmetry(unknown)
+{
+       debugmsg("indexed constructor from ex,ex,ex,ex,ex", LOGLEVEL_CONSTRUCT);
+       tinfo_key = TINFO_indexed;
+       GINAC_ASSERT(all_indices_of_type_idx());
+}
+
+indexed::indexed(const ex & b, symmetry_type symm, const ex & i1, const ex & i2) : inherited(b, i1, i2), symmetry(symm)
 {
-       debugmsg("indexed constructor from ex,ex,ex,ex",LOGLEVEL_CONSTRUCT);
-       tinfo_key=TINFO_indexed;
-       GINAC_ASSERT(all_of_type_idx());
+       debugmsg("indexed constructor from ex,symmetry,ex,ex", LOGLEVEL_CONSTRUCT);
+       tinfo_key = TINFO_indexed;
+       GINAC_ASSERT(all_indices_of_type_idx());
 }
 
-/** Construct indexed object with a specified vector of indices. The indices
- *  must be of class idx or a subclass.
- *
- *  @param iv Vector of indices
- *  @return newly constructed indexed object */
-indexed::indexed(const exvector & iv) : inherited(iv)
+indexed::indexed(const ex & b, symmetry_type symm, const ex & i1, const ex & i2, const ex & i3) : inherited(b, i1, i2, i3), symmetry(symm)
 {
-       debugmsg("indexed constructor from exvector",LOGLEVEL_CONSTRUCT);
-       tinfo_key=TINFO_indexed;
-       GINAC_ASSERT(all_of_type_idx());
+       debugmsg("indexed constructor from ex,symmetry,ex,ex,ex", LOGLEVEL_CONSTRUCT);
+       tinfo_key = TINFO_indexed;
+       GINAC_ASSERT(all_indices_of_type_idx());
+}
+
+indexed::indexed(const ex & b, symmetry_type symm, const ex & i1, const ex & i2, const ex & i3, const ex & i4) : inherited(b, i1, i2, i3, i4), symmetry(symm)
+{
+       debugmsg("indexed constructor from ex,symmetry,ex,ex,ex,ex", LOGLEVEL_CONSTRUCT);
+       tinfo_key = TINFO_indexed;
+       GINAC_ASSERT(all_indices_of_type_idx());
+}
+
+indexed::indexed(const ex & b, const exvector & v) : inherited(b), symmetry(unknown)
+{
+       debugmsg("indexed constructor from ex,exvector", LOGLEVEL_CONSTRUCT);
+       seq.insert(seq.end(), v.begin(), v.end());
+       tinfo_key = TINFO_indexed;
+       GINAC_ASSERT(all_indices_of_type_idx());
+}
+
+indexed::indexed(const ex & b, symmetry_type symm, const exvector & v) : inherited(b), symmetry(symm)
+{
+       debugmsg("indexed constructor from ex,symmetry,exvector", LOGLEVEL_CONSTRUCT);
+       seq.insert(seq.end(), v.begin(), v.end());
+       tinfo_key = TINFO_indexed;
+       GINAC_ASSERT(all_indices_of_type_idx());
+}
+
+indexed::indexed(symmetry_type symm, const exprseq & es) : inherited(es), symmetry(symm)
+{
+       debugmsg("indexed constructor from symmetry,exprseq", LOGLEVEL_CONSTRUCT);
+       tinfo_key = TINFO_indexed;
+       GINAC_ASSERT(all_indices_of_type_idx());
 }
 
-indexed::indexed(exvector * ivp) : inherited(ivp)
+indexed::indexed(symmetry_type symm, const exvector & v, bool discardable) : inherited(v, discardable), symmetry(symm)
 {
-       debugmsg("indexed constructor from exvector *",LOGLEVEL_CONSTRUCT);
-       tinfo_key=TINFO_indexed;
-       GINAC_ASSERT(all_of_type_idx());
+       debugmsg("indexed constructor from symmetry,exvector", LOGLEVEL_CONSTRUCT);
+       tinfo_key = TINFO_indexed;
+       GINAC_ASSERT(all_indices_of_type_idx());
+}
+
+indexed::indexed(symmetry_type symm, exvector * vp) : inherited(vp), symmetry(symm)
+{
+       debugmsg("indexed constructor from symmetry,exvector *", LOGLEVEL_CONSTRUCT);
+       tinfo_key = TINFO_indexed;
+       GINAC_ASSERT(all_indices_of_type_idx());
 }
 
 //////////
@@ -146,7 +163,9 @@ indexed::indexed(exvector * ivp) : inherited(ivp)
 indexed::indexed(const archive_node &n, const lst &sym_lst) : inherited(n, sym_lst)
 {
        debugmsg("indexed constructor from archive_node", LOGLEVEL_CONSTRUCT);
-       tinfo_key = TINFO_indexed;
+       unsigned int symm;
+       if (!(n.find_unsigned("symmetry", symm)))
+               throw (std::runtime_error("unknown indexed symmetry type in archive"));
 }
 
 /** Unarchive the object. */
@@ -159,90 +178,202 @@ ex indexed::unarchive(const archive_node &n, const lst &sym_lst)
 void indexed::archive(archive_node &n) const
 {
        inherited::archive(n);
+       n.add_unsigned("symmetry", symmetry);
 }
 
 //////////
 // functions overriding virtual functions from bases classes
 //////////
 
-// public
-
 void indexed::printraw(std::ostream & os) const
 {
-       debugmsg("indexed printraw",LOGLEVEL_PRINT);
-       os << "indexed(indices=";
+       debugmsg("indexed printraw", 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);
-       os << std::string(indent,' ') << "indexed: " << seq.size() << " indices";
+       debugmsg("indexed printtree", LOGLEVEL_PRINT);
+       GINAC_ASSERT(seq.size() > 0);
+
+       os << std::string(indent, ' ') << class_name() << ", " << seq.size()-1 << " indices";
        os << ",hash=" << hashvalue << ",flags=" << flags << std::endl;
-       printtreeindices(os,indent);
+       printtreeindices(os, indent);
 }
 
 void indexed::print(std::ostream & os, unsigned upper_precedence) const
 {
-       debugmsg("indexed print",LOGLEVEL_PRINT);
-       os << "UNNAMEDINDEX";
+       debugmsg("indexed print", LOGLEVEL_PRINT);
+       GINAC_ASSERT(seq.size() > 0);
+
+       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);
 }
 
-void indexed::printcsrc(std::ostream & os, unsigned type,
-                        unsigned upper_precedence) const
-{
-       debugmsg("indexed print csrc",LOGLEVEL_PRINT);
-       print(os,upper_precedence);
-}
-
 bool indexed::info(unsigned inf) const
 {
-       if (inf==info_flags::indexed) return true;
-       if (inf==info_flags::has_indices) return seq.size()!=0;
+       if (inf == info_flags::indexed) return true;
+       if (inf == info_flags::has_indices) return seq.size() > 1;
        return inherited::info(inf);
 }
 
-// protected
-
-/** Implementation of ex::diff() for an indexed object. It always returns 0.
- *  @see ex::diff */
-ex indexed::derivative(const symbol & s) const
+bool indexed::all_index_values_are(unsigned inf) const
 {
-       return _ex0();
+       // No indices? Then no property can be fulfilled
+       if (seq.size() < 2)
+               return false;
+
+       // Check all indices
+       exvector::const_iterator it = seq.begin() + 1, itend = seq.end();
+       while (it != itend) {
+               GINAC_ASSERT(is_ex_of_type(*it, idx));
+               if (!ex_to_idx(*it).get_value().info(inf))
+                       return false;
+               it++;
+       }
+       return true;
 }
 
 int indexed::compare_same_type(const basic & other) const
 {
-       GINAC_ASSERT(is_of_type(other,indexed));
+       GINAC_ASSERT(is_of_type(other, indexed));
        return inherited::compare_same_type(other);
 }
 
-bool indexed::is_equal_same_type(const basic & other) const
+// The main difference between sort_index_vector() and canonicalize_indices()
+// is that the latter takes the symmetry of the object into account. Once we
+// implement mixed symmetries, canonicalize_indices() will only be able to
+// reorder index pairs with known symmetry properties, while sort_index_vector()
+// always sorts the whole vector.
+
+/** Bring a vector of indices into a canonic order (don't care about the
+ *  symmetry of the objects carrying the indices). Dummy indices will lie
+ *  next to each other after the sorting.
+ *
+ *  @param v Index vector to be sorted */
+static void sort_index_vector(exvector &v)
 {
-       GINAC_ASSERT(is_of_type(other,indexed));
-       return inherited::is_equal_same_type(other);
+       // Nothing to sort if less than 2 elements
+       if (v.size() < 2)
+               return;
+
+       // Simple bubble sort algorithm should be sufficient for the small
+       // number of indices expected
+       exvector::iterator it1 = v.begin(), itend = v.end(), next_to_last_idx = itend - 1;
+       while (it1 != next_to_last_idx) {
+               exvector::iterator it2 = it1 + 1;
+               while (it2 != itend) {
+                       if (it1->compare(*it2) > 0)
+                               it1->swap(*it2);
+                       it2++;
+               }
+               it1++;
+       }
 }
 
-unsigned indexed::return_type(void) const
+/** Bring a vector of indices into a canonic order. This operation only makes
+ *  sense if the object carrying these indices is either symmetric or totally
+ *  antisymmetric with respect to the indices.
+ *
+ *  @param itbegin Start of index vector
+ *  @param itend End of index vector
+ *  @param antisymm Whether the object is antisymmetric
+ *  @return the sign introduced by the reordering of the indices if the object
+ *          is antisymmetric (or 0 if two equal indices are encountered). For
+ *          symmetric objects, this is always +1. If the index vector was
+ *          already in a canonic order this function returns INT_MAX. */
+static int canonicalize_indices(exvector::iterator itbegin, exvector::iterator itend, bool antisymm)
 {
-       return return_types::noncommutative;
+       bool something_changed = false;
+       int sig = 1;
+
+       // Simple bubble sort algorithm should be sufficient for the small
+       // number of indices expected
+       exvector::iterator it1 = itbegin, next_to_last_idx = itend - 1;
+       while (it1 != next_to_last_idx) {
+               exvector::iterator it2 = it1 + 1;
+               while (it2 != itend) {
+                       int cmpval = it1->compare(*it2);
+                       if (cmpval == 1) {
+                               it1->swap(*it2);
+                               something_changed = true;
+                               if (antisymm)
+                                       sig = -sig;
+                       } else if (cmpval == 0 && antisymm) {
+                               something_changed = true;
+                               sig = 0;
+                       }
+                       it2++;
+               }
+               it1++;
+       }
+
+       return something_changed ? sig : INT_MAX;
 }
-   
-unsigned indexed::return_type_tinfo(void) const
+
+ex indexed::eval(int level) const
 {
-       return tinfo_key;
+       // First evaluate children, then we will end up here again
+       if (level > 1)
+               return indexed(symmetry, evalchildren(level));
+
+       // Canonicalize indices according to the symmetry properties
+       if (seq.size() > 2 && (symmetry != unknown && symmetry != mixed)) {
+               exvector v = seq;
+               int sig = canonicalize_indices(v.begin() + 1, v.end(), symmetry == antisymmetric);
+               if (sig != INT_MAX) {
+                       // Something has changed while sorting indices, more evaluations later
+                       if (sig == 0)
+                               return _ex0();
+                       return ex(sig) * thisexprseq(v);
+               }
+       }
+
+       // Let the class of the base object perform additional evaluations
+       return op(0).bp->eval_indexed(*this);
 }
 
 ex indexed::thisexprseq(const exvector & v) const
 {
-       return indexed(v);
+       return indexed(symmetry, v);
 }
 
 ex indexed::thisexprseq(exvector * vp) const
 {
-       return indexed(vp);
+       return indexed(symmetry, vp);
+}
+
+ex indexed::expand(unsigned options) const
+{
+       GINAC_ASSERT(seq.size() > 0);
+
+       if ((options & expand_options::expand_indexed) && is_ex_exactly_of_type(seq[0], add)) {
+
+               // expand_indexed expands (a+b).i -> a.i + b.i
+               const ex & base = seq[0];
+               ex sum = _ex0();
+               for (unsigned i=0; i<base.nops(); i++) {
+                       exvector s = seq;
+                       s[0] = base.op(i);
+                       sum += thisexprseq(s).expand();
+               }
+               return sum;
+
+       } else
+               return inherited::expand(options);
 }
 
 //////////
@@ -255,57 +386,391 @@ ex indexed::thisexprseq(exvector * vp) const
 // non-virtual functions in this class
 //////////
 
-// protected
-
 void indexed::printrawindices(std::ostream & os) const
 {
-       if (seq.size()!=0) {
-               for (exvector::const_iterator cit=seq.begin(); cit!=seq.end(); ++cit) {
-                       (*cit).printraw(os);
-                       os << ",";
+       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()!=0) {
-               for (exvector::const_iterator cit=seq.begin(); cit!=seq.end(); ++cit) {
-                       os << std::string(indent+delta_indent,' ');
-                       (*cit).printraw(os);
+       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
 {
-       if (seq.size()!=0) {
-               if (seq.size()>1) {
-                       os << "{";
-               }
-               exvector::const_iterator last=seq.end()-1;
-               exvector::const_iterator cit=seq.begin();
-               for (; cit!=last; ++cit) {
-                       (*cit).print(os);
-                       os << ",";
+       if (seq.size() > 1) {
+               exvector::const_iterator it=seq.begin() + 1, itend = seq.end();
+               while (it != itend) {
+                       it->print(os);
+                       it++;
                }
-               (*cit).print(os);
-               if (seq.size()>1) {
-                       os << "}";
+       }
+}
+
+/** Check whether all indices are of class idx. This function is used
+ *  internally to make sure that all constructed indexed objects really
+ *  carry indices and not some other classes. */
+bool indexed::all_indices_of_type_idx(void) const
+{
+       GINAC_ASSERT(seq.size() > 0);
+       exvector::const_iterator it = seq.begin() + 1, itend = seq.end();
+       while (it != itend) {
+               if (!is_ex_of_type(*it, idx))
+                       return false;
+               it++;
+       }
+       return true;
+}
+
+//////////
+// 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 all indices are of class idx or a subclass. This function
- *  is used internally to make sure that all constructed indexed objects
- *  really carry indices and not some other classes. */
-bool indexed::all_of_type_idx(void) const
+/** Check whether two sorted index vectors are consistent (i.e. equal). */
+static bool indices_consistent(const exvector & v1, const exvector & v2)
 {
-       for (exvector::const_iterator cit=seq.begin(); cit!=seq.end(); ++cit) {
-               if (!is_ex_of_type(*cit,idx)) return false;
+       // Number of indices must be the same
+       if (v1.size() != v2.size())
+               return false;
+
+       // And also the indices themselves
+       exvector::const_iterator ait = v1.begin(), aitend = v1.end(),
+                                bit = v2.begin(), bitend = v2.end();
+       while (ait != aitend) {
+               if (!ait->is_equal(*bit))
+                       return false;
+               ait++; bit++;
        }
        return true;
 }
 
+exvector indexed::get_free_indices(void) const
+{
+       exvector free_indices, dummy_indices;
+       find_free_and_dummy(seq.begin() + 1, seq.end(), free_indices, dummy_indices);
+       return free_indices;
+}
+
+exvector add::get_free_indices(void) const
+{
+       exvector free_indices;
+       for (unsigned i=0; i<nops(); i++) {
+               if (i == 0)
+                       free_indices = op(i).get_free_indices();
+               else {
+                       exvector free_indices_of_term = op(i).get_free_indices();
+                       if (!indices_consistent(free_indices, free_indices_of_term))
+                               throw (std::runtime_error("add::get_free_indices: inconsistent indices in sum"));
+               }
+       }
+       return free_indices;
+}
+
+exvector mul::get_free_indices(void) const
+{
+       // Concatenate free indices of all factors
+       exvector un;
+       for (unsigned i=0; i<nops(); i++) {
+               exvector free_indices_of_factor = op(i).get_free_indices();
+               un.insert(un.end(), free_indices_of_factor.begin(), free_indices_of_factor.end());
+       }
+
+       // And remove the dummy indices
+       exvector free_indices, dummy_indices;
+       find_free_and_dummy(un.begin(), un.end(), free_indices, dummy_indices);
+       return free_indices;
+}
+
+exvector ncmul::get_free_indices(void) const
+{
+       // Concatenate free indices of all factors
+       exvector un;
+       for (unsigned i=0; i<nops(); i++) {
+               exvector free_indices_of_factor = op(i).get_free_indices();
+               un.insert(un.end(), free_indices_of_factor.begin(), free_indices_of_factor.end());
+       }
+
+       // And remove the dummy indices
+       exvector free_indices, dummy_indices;
+       find_free_and_dummy(un.begin(), un.end(), free_indices, dummy_indices);
+       return free_indices;
+}
+
+exvector power::get_free_indices(void) const
+{
+       // Return free indices of basis
+       return basis.get_free_indices();
+}
+
+/** Simplify product of indexed expressions (commutative, noncommutative and
+ *  simple squares), return list of free indices. */
+ex simplify_indexed_product(const ex & e, exvector & free_indices, const scalar_products & sp)
+{
+       // Remember whether the product was commutative or noncommutative
+       // (because we chop it into factors and need to reassemble later)
+       bool non_commutative = is_ex_exactly_of_type(e, ncmul);
+
+       // Collect factors in an exvector, store squares twice
+       exvector v;
+       v.reserve(e.nops() * 2);
+
+       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()));
+               v.push_back(e.op(0));
+               v.push_back(e.op(0));
+       } else {
+               for (int 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())) {
+                               v.push_back(f.op(0));
+                   v.push_back(f.op(0));
+                       } 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++)
+                                       v.push_back(f.op(j));
+                       } else
+                               v.push_back(f);
+               }
+       }
+
+       // Perform contractions
+       bool something_changed = false;
+       GINAC_ASSERT(v.size() > 1);
+       exvector::iterator it1, itend = v.end(), next_to_last = itend - 1;
+       for (it1 = v.begin(); it1 != next_to_last; it1++) {
+
+try_again:
+               if (!is_ex_of_type(*it1, indexed))
+                       continue;
+
+               // Indexed factor found, look for contraction candidates
+               exvector::iterator it2;
+               for (it2 = it1 + 1; it2 != itend; it2++) {
+
+                       if (!is_ex_of_type(*it2, indexed))
+                               continue;
+
+                       // 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);
+                       if (dummy.size() == 0)
+                               continue;
+
+                       // At least one dummy index, is it a defined scalar product?
+                       if (free.size() == 0) {
+                               if (sp.is_defined(*it1, *it2)) {
+                                       *it1 = sp.evaluate(*it1, *it2);
+                                       *it2 = _ex1();
+                                       something_changed = true;
+                                       goto try_again;
+                               }
+                       }
+
+                       // Try to contract the first one with the second one
+                       bool contracted = it1->op(0).bp->contract_with(*it1, *it2);
+                       if (!contracted) {
+
+                               // That didn't work; maybe the second object knows how to
+                               // contract itself with the first one
+                               contracted = it2->op(0).bp->contract_with(*it2, *it1);
+                       }
+                       if (contracted) {
+                               something_changed = true;
+
+                               // Both objects may have new indices now or they might
+                               // even not be indexed objects any more, so we have to
+                               // start over
+                               goto try_again;
+                       }
+               }
+       }
+
+       // Find free indices (concatenate them all and call find_free_and_dummy())
+       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());
+               }
+               it1++;
+       }
+       find_free_and_dummy(un.begin(), un.end(), free_indices, dummy_indices);
+
+       if (something_changed) {
+               if (non_commutative)
+                       return ncmul(v);
+               else
+                       return mul(v);
+       } else
+               return e;
+}
+
+/** Simplify indexed expression, return list of free indices. */
+ex simplify_indexed(const ex & e, exvector & free_indices, const scalar_products & sp)
+{
+       // Expand the expression
+       ex e_expanded = e.expand();
+
+       // Simplification of single indexed object: just find the free indices
+       if (is_ex_of_type(e_expanded, indexed)) {
+               const indexed &i = ex_to_indexed(e_expanded);
+               exvector dummy_indices;
+               find_free_and_dummy(i.seq.begin() + 1, i.seq.end(), free_indices, dummy_indices);
+               return e_expanded;
+       }
+
+       // Simplification of sum = sum of simplifications, check consistency of
+       // free indices in each term
+       if (is_ex_exactly_of_type(e_expanded, add)) {
+               ex sum = _ex0();
+
+               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"));
+               }
+
+               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, sp);
+
+       // Cannot do anything
+       free_indices.clear();
+       return e_expanded;
+}
+
+ex simplify_indexed(const ex & e)
+{
+       exvector free_indices;
+       scalar_products sp;
+       return simplify_indexed(e, free_indices, sp);
+}
+
+ex simplify_indexed(const ex & e, const scalar_products & sp)
+{
+       exvector free_indices;
+       return simplify_indexed(e, free_indices, sp);
+}
+
+//////////
+// helper classes
+//////////
+
+void scalar_products::add(const ex & v1, const ex & v2, const ex & sp)
+{
+       spm[make_key(v1, v2)] = sp;
+}
+
+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;
+       for (spmap::const_iterator cit=spm.begin(); cit!=spm.end(); ++cit) {
+               const spmapkey & k = cit->first;
+               std::cerr << "item key=(" << k.first << "," << k.second;
+               std::cerr << "), value=" << cit->second << std::endl;
+       }
+}
+
+/** 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);
+}
+
 } // namespace GiNaC
index 94ba0e6..b15ed8b 100644 (file)
@@ -1,6 +1,6 @@
 /** @file indexed.h
  *
- *  Interface to GiNaC's index carrying objects. */
+ *  Interface to GiNaC's indexed expressions. */
 
 /*
  *  GiNaC Copyright (C) 1999-2001 Johannes Gutenberg University Mainz, Germany
 #ifndef __GINAC_INDEXED_H__
 #define __GINAC_INDEXED_H__
 
-#include <string>
+#include <map>
+
 #include "exprseq.h"
 
 namespace GiNaC {
 
 
-/** Base class for objects with indices. */
+class scalar_products;
+
+/** This class holds an indexed expression. It consists of a 'base' expression
+ *  (the expression being indexed) which can be accessed as op(0), and n (n >= 0)
+ *  indices (all of class idx), accessible as op(1)..op(n). */
 class indexed : public exprseq
 {
        GINAC_DECLARE_REGISTERED_CLASS(indexed, exprseq)
 
+       friend ex simplify_indexed(const ex & e, exvector & free_indices, const scalar_products & sp);
+       friend ex simplify_indexed_product(const ex & e, exvector & free_indices, const scalar_products & sp);
+
+       // types
+public:
+       /** Type of symmetry of the object with respect to commutation of its indices. */
+       typedef enum {
+               unknown,       /**< symmetry properties unknown */
+               symmetric,     /**< totally symmetric */
+               antisymmetric, /**< totally antisymmetric */
+               mixed          /**< mixed symmetry (unimplemented) */
+       } symmetry_type;
+
        // other constructors
 public:
-       indexed(const ex & i1);
-       indexed(const ex & i1, const ex & i2);
-       indexed(const ex & i1, const ex & i2, const ex & i3);
-       indexed(const ex & i1, const ex & i2, const ex & i3, const ex & i4);
-       indexed(const exvector & iv);
-       indexed(exvector * iv);
+       /** Construct indexed object with no index.
+        *
+        *  @param b Base expression
+        *  @return newly constructed indexed object */
+       indexed(const ex & b);
+
+       /** Construct indexed object with one index. The index must be of class idx.
+        *
+        *  @param b Base expression
+        *  @param i1 The index
+        *  @return newly constructed indexed object */
+       indexed(const ex & b, const ex & i1);
+
+       /** Construct indexed object with two indices. The indices must be of class idx.
+        *
+        *  @param b Base expression
+        *  @param i1 First index
+        *  @param i2 Second index
+        *  @return newly constructed indexed object */
+       indexed(const ex & b, const ex & i1, const ex & i2);
+
+       /** Construct indexed object with three indices. The indices must be of class idx.
+        *
+        *  @param b Base expression
+        *  @param i1 First index
+        *  @param i2 Second index
+        *  @param i3 Third index
+        *  @return newly constructed indexed object */
+       indexed(const ex & b, const ex & i1, const ex & i2, const ex & i3);
+
+       /** Construct indexed object with four indices. The indices must be of class idx.
+        *
+        *  @param b Base expression
+        *  @param i1 First index
+        *  @param i2 Second index
+        *  @param i3 Third index
+        *  @param i4 Fourth index
+        *  @return newly constructed indexed object */
+       indexed(const ex & b, const ex & i1, const ex & i2, const ex & i3, const ex & i4);
+
+       /** Construct indexed object with two indices and a specified symmetry. The
+        *  indices must be of class idx.
+        *
+        *  @param b Base expression
+        *  @param symm Symmetry of indices
+        *  @param i1 First index
+        *  @param i2 Second index
+        *  @return newly constructed indexed object */
+       indexed(const ex & b, symmetry_type symm, const ex & i1, const ex & i2);
+
+       /** Construct indexed object with three indices and a specified symmetry.
+        *  The indices must be of class idx.
+        *
+        *  @param b Base expression
+        *  @param symm Symmetry of indices
+        *  @param i1 First index
+        *  @param i2 Second index
+        *  @param i3 Third index
+        *  @return newly constructed indexed object */
+       indexed(const ex & b, symmetry_type symm, const ex & i1, const ex & i2, const ex & i3);
+
+       /** Construct indexed object with four indices and a specified symmetry. The
+        *  indices must be of class idx.
+        *
+        *  @param b Base expression
+        *  @param symm Symmetry of indices
+        *  @param i1 First index
+        *  @param i2 Second index
+        *  @param i3 Third index
+        *  @param i4 Fourth index
+        *  @return newly constructed indexed object */
+       indexed(const ex & b, symmetry_type symm, const ex & i1, const ex & i2, const ex & i3, const ex & i4);
+
+       /** Construct indexed object with a specified vector of indices. The indices
+        *  must be of class idx.
+        *
+        *  @param b Base expression
+        *  @param iv Vector of indices
+        *  @return newly constructed indexed object */
+       indexed(const ex & b, const exvector & iv);
+
+       /** Construct indexed object with a specified vector of indices and
+        *  symmetry. The indices must be of class idx.
+        *
+        *  @param b Base expression
+        *  @param symm Symmetry of indices
+        *  @param iv Vector of indices
+        *  @return newly constructed indexed object */
+       indexed(const ex & b, symmetry_type symm, const exvector & iv);
+
+       // internal constructors
+       indexed(symmetry_type symm, const exprseq & es);
+       indexed(symmetry_type symm, const exvector & v, bool discardable = false);
+       indexed(symmetry_type symm, exvector * vp); // vp will be deleted
 
        // functions overriding virtual functions from base classes
 public:
        void printraw(std::ostream & os) const;
        void printtree(std::ostream & os, unsigned indent) const;
        void print(std::ostream & os, unsigned upper_precedence=0) const;
-       void printcsrc(std::ostream & os, unsigned type, unsigned upper_precedence) const;
        bool info(unsigned inf) const;
-
-       /** Return the vector of indices on this object. */
-       exvector get_indices(void) const {return seq;}
+       ex eval(int level = 0) const;
+       exvector get_free_indices(void) const;
 
 protected:
-       ex derivative(const symbol & s) const;
-       bool is_equal_same_type(const basic & other) const;
-       unsigned return_type(void) const;
-       unsigned return_type_tinfo(void) const;
        ex thisexprseq(const exvector & v) const;
        ex thisexprseq(exvector * vp) const;
+       unsigned return_type(void) const { return return_types::commutative; }
+       ex expand(unsigned options = 0) const;
 
        // new virtual functions which can be overridden by derived classes
        // none
        
        // non-virtual functions in this class
+public:
+       /** Check whether all index values have a certain property.
+        *  @see class info_flags */
+       bool all_index_values_are(unsigned inf) const;
+
 protected:
        void printrawindices(std::ostream & os) const;
        void printtreeindices(std::ostream & os, unsigned indent) const;
        void printindices(std::ostream & os) const;
-       bool all_of_type_idx(void) const;
+       bool all_indices_of_type_idx(void) const;
 
-// member variables
-       // none
+       // member variables
+protected:
+       symmetry_type symmetry; /**< Index symmetry */
 };
 
+
+typedef std::pair<ex, ex> spmapkey;
+
+struct spmapkey_is_less {
+       bool operator() (const spmapkey &p, const spmapkey &q) const 
+       {
+               int cmp = p.first.compare(q.first);
+               return ((cmp<0) || (!(cmp>0) && p.second.compare(q.second)<0));
+       }
+};
+
+typedef std::map<spmapkey, ex, spmapkey_is_less> spmap;
+
+/** Helper class for storing information about known scalar products which
+ *  are to be automatically replaced by simplify_indexed().
+ *
+ *  @see simplify_indexed */
+class scalar_products {
+public:
+       /** Register scalar product pair and its value. */
+       void add(const ex & v1, const ex & v2, const ex & sp);
+
+       /** Clear all registered scalar products. */
+       void clear(void);
+
+       bool is_defined(const ex & v1, const ex & v2) const;
+       ex evaluate(const ex & v1, const ex & v2) const;
+       void debugprint(void) const;
+
+private:
+       static spmapkey make_key(const ex & v1, const ex & v2);
+
+       spmap spm; /*< Map from defined scalar product pairs to their values */
+};
+
+
 // utility functions
 inline const indexed &ex_to_indexed(const ex &e)
 {
        return static_cast<const indexed &>(*e.bp);
 }
 
+
+/** Simplify/canonicalize expression containing indexed objects. This
+ *  performs contraction of dummy indices where possible and checks whether
+ *  the free indices in sums are consistent.
+ *
+ *  @param e The expression to be simplified
+ *  @return simplified expression */
+ex simplify_indexed(const ex & e);
+
+/** 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 e The expression to be simplified
+ *  @param sp Scalar products to be replaced automatically
+ *  @return simplified expression */
+ex simplify_indexed(const ex & e, const scalar_products & sp);
+
+
 } // namespace GiNaC
 
 #endif // ndef __GINAC_INDEXED_H__
index 0c32a76..a63275f 100644 (file)
@@ -28,6 +28,8 @@
 #include "archive.h"
 #include "numeric.h"
 #include "lst.h"
+#include "idx.h"
+#include "indexed.h"
 #include "utils.h"
 #include "debugmsg.h"
 #include "power.h"
@@ -157,7 +159,7 @@ void matrix::print(std::ostream & os, unsigned upper_precedence) const
 void matrix::printraw(std::ostream & os) const
 {
        debugmsg("matrix printraw",LOGLEVEL_PRINT);
-       os << "matrix(" << row << "," << col <<",";
+       os << class_name() << "(" << row << "," << col <<",";
        for (unsigned r=0; r<row-1; ++r) {
                os << "(";
                for (unsigned c=0; c<col-1; ++c)
@@ -217,7 +219,7 @@ bool matrix::has(const ex & other) const
        return false;
 }
 
-/** evaluate matrix entry by entry. */
+/** Evaluate matrix entry by entry. */
 ex matrix::eval(int level) const
 {
        debugmsg("matrix eval",LOGLEVEL_MEMBER_FUNCTION);
@@ -241,7 +243,7 @@ ex matrix::eval(int level) const
                                                                                           status_flags::evaluated );
 }
 
-/** evaluate matrix numerically entry by entry. */
+/** Evaluate matrix numerically entry by entry. */
 ex matrix::evalf(int level) const
 {
        debugmsg("matrix evalf",LOGLEVEL_MEMBER_FUNCTION);
@@ -265,6 +267,16 @@ ex matrix::evalf(int level) const
        return matrix(row, col, m2);
 }
 
+ex matrix::subs(const lst & ls, const lst & lr) const
+{
+       exvector m2(row * col);
+       for (unsigned r=0; r<row; ++r)
+               for (unsigned c=0; c<col; ++c)
+                       m2[r*col+c] = m[r*col+c].subs(ls, lr);
+
+       return matrix(row, col, m2);
+}
+
 // protected
 
 int matrix::compare_same_type(const basic & other) const
@@ -292,6 +304,189 @@ int matrix::compare_same_type(const basic & other) const
        return 0;
 }
 
+/** Automatic symbolic evaluation of an indexed matrix. */
+ex matrix::eval_indexed(const basic & i) const
+{
+       GINAC_ASSERT(is_of_type(i, indexed));
+       GINAC_ASSERT(is_ex_of_type(i.op(0), matrix));
+
+       bool all_indices_unsigned = static_cast<const indexed &>(i).all_index_values_are(info_flags::nonnegint);
+
+       // Check indices
+       if (i.nops() == 2) {
+
+               // One index, must be one-dimensional vector
+               if (row != 1 && col != 1)
+                       throw (std::runtime_error("matrix::eval_indexed(): vector must have exactly 1 index"));
+
+               const idx & i1 = ex_to_idx(i.op(1));
+
+               if (col == 1) {
+
+                       // Column vector
+                       if (!i1.get_dim().is_equal(row))
+                               throw (std::runtime_error("matrix::eval_indexed(): dimension of index must match number of vector elements"));
+
+                       // Index numeric -> return vector element
+                       if (all_indices_unsigned) {
+                               unsigned n1 = ex_to_numeric(i1.get_value()).to_int();
+                               if (n1 >= row)
+                                       throw (std::runtime_error("matrix::eval_indexed(): value of index exceeds number of vector elements"));
+                               return (*this)(n1, 0);
+                       }
+
+               } else {
+
+                       // Row vector
+                       if (!i1.get_dim().is_equal(col))
+                               throw (std::runtime_error("matrix::eval_indexed(): dimension of index must match number of vector elements"));
+
+                       // Index numeric -> return vector element
+                       if (all_indices_unsigned) {
+                               unsigned n1 = ex_to_numeric(i1.get_value()).to_int();
+                               if (n1 >= col)
+                                       throw (std::runtime_error("matrix::eval_indexed(): value of index exceeds number of vector elements"));
+                               return (*this)(0, n1);
+                       }
+               }
+
+       } else if (i.nops() == 3) {
+
+               // Two indices
+               const idx & i1 = ex_to_idx(i.op(1));
+               const idx & i2 = ex_to_idx(i.op(2));
+
+               if (!i1.get_dim().is_equal(row))
+                       throw (std::runtime_error("matrix::eval_indexed(): dimension of first index must match number of rows"));
+               if (!i2.get_dim().is_equal(col))
+                       throw (std::runtime_error("matrix::eval_indexed(): dimension of second index must match number of columns"));
+
+               // Pair of dummy indices -> compute trace
+               if (is_dummy_pair(i1, i2))
+                       return trace();
+
+               // Both indices numeric -> return matrix element
+               if (all_indices_unsigned) {
+                       unsigned n1 = ex_to_numeric(i1.get_value()).to_int(), n2 = ex_to_numeric(i2.get_value()).to_int();
+                       if (n1 >= row)
+                               throw (std::runtime_error("matrix::eval_indexed(): value of first index exceeds number of rows"));
+                       if (n2 >= col)
+                               throw (std::runtime_error("matrix::eval_indexed(): value of second index exceeds number of columns"));
+                       return (*this)(n1, n2);
+               }
+
+       } else
+               throw (std::runtime_error("matrix::eval_indexed(): matrix must have exactly 2 indices"));
+
+       return i.hold();
+}
+
+/** Contraction of an indexed matrix with something else. */
+bool matrix::contract_with(ex & self, ex & other) const
+{
+       GINAC_ASSERT(is_ex_of_type(self, indexed));
+       GINAC_ASSERT(is_ex_of_type(other, indexed));
+       GINAC_ASSERT(self.nops() == 2 || self.nops() == 3);
+       GINAC_ASSERT(is_ex_of_type(self.op(0), matrix));
+
+       // Only contract with other matrices
+       if (!is_ex_of_type(other.op(0), matrix))
+               return false;
+
+       GINAC_ASSERT(other.nops() == 2 || other.nops() == 3);
+
+       const matrix &self_matrix = ex_to_matrix(self.op(0));
+       const matrix &other_matrix = ex_to_matrix(other.op(0));
+
+       if (self.nops() == 2) {
+               unsigned self_dim = (self_matrix.col == 1) ? self_matrix.row : self_matrix.col;
+
+               if (other.nops() == 2) { // vector * vector (scalar product)
+                       unsigned other_dim = (other_matrix.col == 1) ? other_matrix.row : other_matrix.col;
+
+                       if (self_matrix.col == 1) {
+                               if (other_matrix.col == 1) {
+                                       // Column vector * column vector, transpose first vector
+                                       self = self_matrix.transpose().mul(other_matrix)(0, 0);
+                               } else {
+                                       // Column vector * row vector, swap factors
+                                       self = other_matrix.mul(self_matrix)(0, 0);
+                               }
+                       } else {
+                               if (other_matrix.col == 1) {
+                                       // Row vector * column vector, perfect
+                                       self = self_matrix.mul(other_matrix)(0, 0);
+                               } else {
+                                       // Row vector * row vector, transpose second vector
+                                       self = self_matrix.mul(other_matrix.transpose())(0, 0);
+                               }
+                       }
+                       other = _ex1();
+                       return true;
+
+               } else { // vector * matrix
+
+                       GINAC_ASSERT(other.nops() == 3);
+
+                       // B_i * A_ij = (B*A)_j (B is row vector)
+                       if (is_dummy_pair(self.op(1), other.op(1))) {
+                               if (self_matrix.row == 1)
+                                       self = indexed(self_matrix.mul(other_matrix), other.op(2));
+                               else
+                                       self = indexed(self_matrix.transpose().mul(other_matrix), other.op(2));
+                               other = _ex1();
+                               return true;
+                       }
+
+                       // B_j * A_ij = (A*B)_i (B is column vector)
+                       if (is_dummy_pair(self.op(1), other.op(2))) {
+                               if (self_matrix.col == 1)
+                                       self = indexed(other_matrix.mul(self_matrix), other.op(1));
+                               else
+                                       self = indexed(other_matrix.mul(self_matrix.transpose()), other.op(1));
+                               other = _ex1();
+                               return true;
+                       }
+               }
+
+       } else if (other.nops() == 3) { // matrix * matrix
+
+               GINAC_ASSERT(self.nops() == 3);
+               GINAC_ASSERT(other.nops() == 3);
+
+               // A_ij * B_jk = (A*B)_ik
+               if (is_dummy_pair(self.op(2), other.op(1))) {
+                       self = indexed(self_matrix.mul(other_matrix), self.op(1), other.op(2));
+                       other = _ex1();
+                       return true;
+               }
+
+               // A_ij * B_kj = (A*Btrans)_ik
+               if (is_dummy_pair(self.op(2), other.op(2))) {
+                       self = indexed(self_matrix.mul(other_matrix.transpose()), self.op(1), other.op(1));
+                       other = _ex1();
+                       return true;
+               }
+
+               // A_ji * B_jk = (Atrans*B)_ik
+               if (is_dummy_pair(self.op(1), other.op(1))) {
+                       self = indexed(self_matrix.transpose().mul(other_matrix), self.op(2), other.op(2));
+                       other = _ex1();
+                       return true;
+               }
+
+               // A_ji * B_kj = (B*A)_ki
+               if (is_dummy_pair(self.op(1), other.op(2))) {
+                       self = indexed(other_matrix.mul(self_matrix), other.op(1), self.op(2));
+                       other = _ex1();
+                       return true;
+               }
+       }
+
+       return false;
+}
+
+
 //////////
 // non-virtual functions in this class
 //////////
index 02203e0..7efbfed 100644 (file)
@@ -50,7 +50,9 @@ public:
        bool has(const ex & other) const;
        ex eval(int level=0) const;
        ex evalf(int level=0) const;
-       // ex subs(const lst & ls, const lst & lr) const;
+       ex subs(const lst & ls, const lst & lr) const;
+       ex eval_indexed(const basic & i) const;
+       bool contract_with(ex & self, ex & other) const;
 protected:
        unsigned return_type(void) const { return return_types::noncommutative; };
        // new virtual functions which can be overridden by derived classes
index 887a916..de5d03c 100644 (file)
@@ -418,19 +418,6 @@ ex mul::evalf(int level) const
        return mul(s,overall_coeff.evalf(level));
 }
 
-exvector mul::get_indices(void) const
-{
-       // return union of indices of factors
-       exvector iv;
-       for (epvector::const_iterator cit=seq.begin(); cit!=seq.end(); ++cit) {
-               exvector subiv=(*cit).rest.get_indices();
-               iv.reserve(iv.size()+subiv.size());
-               for (exvector::const_iterator cit2=subiv.begin(); cit2!=subiv.end(); ++cit2)
-                       iv.push_back(*cit2);
-       }
-       return iv;
-}
-
 ex mul::simplify_ncmul(const exvector & v) const
 {
        throw(std::logic_error("mul::simplify_ncmul() should never have been called!"));
index dba5daf..c531fd6 100644 (file)
@@ -61,7 +61,7 @@ public:
        numeric integer_content(void) const;
        ex smod(const numeric &xi) const;
        numeric max_coefficient(void) const;
-       exvector get_indices(void) const;
+       exvector get_free_indices(void) const;
        ex simplify_ncmul(const exvector & v) const;
 protected:
        ex derivative(const symbol & s) const;
index 745cf96..1433a3a 100644 (file)
@@ -471,20 +471,6 @@ ex ncmul::eval(int level) const
                                                                                  status_flags::evaluated);
 }
 
-exvector ncmul::get_indices(void) const
-{
-       // return union of indices of factors
-       exvector iv;
-       for (exvector::const_iterator cit=seq.begin(); cit!=seq.end(); ++cit) {
-               exvector subiv=(*cit).get_indices();
-               iv.reserve(iv.size()+subiv.size());
-               for (exvector::const_iterator cit2=subiv.begin(); cit2!=subiv.end(); ++cit2) {
-                       iv.push_back(*cit2);
-               }
-       }
-       return iv;
-}
-
 ex ncmul::subs(const lst & ls, const lst & lr) const
 {
        return ncmul(subschildren(ls, lr));
index cc186a0..dd96a89 100644 (file)
@@ -63,7 +63,7 @@ public:
        ex coeff(const symbol & s, int n=1) const;
        ex eval(int level=0) const;
        ex subs(const lst & ls, const lst & lr) const;
-       exvector get_indices(void) const;
+       exvector get_free_indices(void) const;
        ex thisexprseq(const exvector & v) const;
        ex thisexprseq(exvector * vp) const;
 protected:
index b40e71f..3bd24a6 100644 (file)
@@ -449,7 +449,7 @@ void numeric::printraw(std::ostream &os) const
        // The method printraw doesn't do much, it simply uses CLN's operator<<()
        // for output, which is ugly but reliable. e.g: 2+2i
        debugmsg("numeric printraw", LOGLEVEL_PRINT);
-       os << "numeric(" << cln::the<cln::cl_N>(value) << ")";
+       os << class_name() << "(" << cln::the<cln::cl_N>(value) << ")";
 }
 
 
index 5ebeae6..b0f9cd2 100644 (file)
@@ -135,7 +135,7 @@ void power::printraw(std::ostream & os) const
 {
        debugmsg("power printraw",LOGLEVEL_PRINT);
 
-       os << "power(";
+       os << class_name() << "(";
        basis.printraw(os);
        os << ",";
        exponent.printraw(os);
@@ -146,8 +146,8 @@ void power::printtree(std::ostream & os, unsigned indent) const
 {
        debugmsg("power printtree",LOGLEVEL_PRINT);
 
-       os << std::string(indent,' ') << "power: "
-          << "hash=" << hashvalue
+       os << std::string(indent,' ') << class_name()
+          << "hash=" << hashvalue
           << " (0x" << std::hex << hashvalue << std::dec << ")"
           << ", flags=" << flags << std::endl;
        basis.printtree(os, indent+delta_indent);
index 8a9d65f..2c74861 100644 (file)
@@ -64,6 +64,7 @@ public:
        ex subs(const lst & ls, const lst & lr) const;
        ex normal(lst &sym_lst, lst &repl_lst, int level = 0) const;
        ex to_rational(lst &repl_lst) const;
+       exvector get_free_indices(void) const;
        ex simplify_ncmul(const exvector & v) const;
 protected:
        ex derivative(const symbol & s) const;
index 692efbc..3b2ac70 100644 (file)
@@ -175,7 +175,7 @@ void pseries::print(std::ostream &os, unsigned upper_precedence) const
 void pseries::printraw(std::ostream &os) const
 {
        debugmsg("pseries printraw", LOGLEVEL_PRINT);
-       os << "pseries(" << var << ";" << point << ";";
+       os << class_name() << "(" << var << ";" << point << ";";
        for (epvector::const_iterator i=seq.begin(); i!=seq.end(); ++i)
                os << "(" << (*i).rest << "," << (*i).coeff << "),";
        os << ")";
@@ -185,7 +185,7 @@ void pseries::printraw(std::ostream &os) const
 void pseries::printtree(std::ostream & os, unsigned indent) const
 {
        debugmsg("pseries printtree",LOGLEVEL_PRINT);
-       os << std::string(indent,' ') << "pseries " 
+       os << std::string(indent,' ') << class_name()
           << ", hash=" << hashvalue
           << " (0x" << std::hex << hashvalue << std::dec << ")"
           << ", flags=" << flags << std::endl;
@@ -227,6 +227,7 @@ int pseries::compare_same_type(const basic & other) const
                ++it;
                ++o_it;
        }
+
        // so they are equal.
        return 0;
 }
index b22a638..20ff40b 100644 (file)
@@ -143,7 +143,7 @@ void relational::print(std::ostream & os, unsigned upper_precedence) const
 void relational::printraw(std::ostream & os) const
 {
        debugmsg("relational printraw",LOGLEVEL_PRINT);
-       os << "RELATIONAL(";
+       os << class_name() << "(";
        lh.printraw(os);
        os << ",";
        rh.printraw(os);
index 35b53e3..e9c3b02 100644 (file)
@@ -92,14 +92,14 @@ void structure::printraw(std::ostream & os) const
 {
        debugmsg("structure printraw",LOGLEVEL_PRINT);
 
-       os << "structure(hash=" << hashvalue << ",flags=" << flags << ")";
+       os << class_name() << "(hash=" << hashvalue << ",flags=" << flags << ")";
 }
 
 void structure::print(std::ostream & os, unsigned upper_precedence) const
 {
        debugmsg("structure print",LOGLEVEL_PRINT);
 
-       os << "structure()";
+       os << class_name() << "()";
 }
 
 void structure::printtree(std::ostream & os, unsigned indent) const
index e46bf4b..9a734c9 100755 (executable)
@@ -331,13 +331,13 @@ basic * ${STRUCTURE}::duplicate() const
 void ${STRUCTURE}::printraw(ostream & os) const
 {
        debugmsg("${STRUCTURE} printraw",LOGLEVEL_PRINT);
-       os << "${STRUCTURE}()";
+       os << class_name() << "()";
 }
 
 void ${STRUCTURE}::print(ostream & os, unsigned upper_precedence) const
 {
        debugmsg("${STRUCTURE} print",LOGLEVEL_PRINT);
-       os << "${STRUCTURE}()";
+       os << class_name() << "()";
 }
 
 void ${STRUCTURE}::printtree(ostream & os, unsigned indent) const
index 169ee63..4328d19 100644 (file)
@@ -25,7 +25,6 @@
 
 #include "symbol.h"
 #include "lst.h"
-#include "idx.h"
 #include "archive.h"
 #include "debugmsg.h"
 #include "utils.h"
@@ -141,7 +140,7 @@ void symbol::print(std::ostream & os, unsigned upper_precedence) const
 void symbol::printraw(std::ostream & os) const
 {
        debugmsg("symbol printraw",LOGLEVEL_PRINT);
-       os << "symbol(" << "name=" << name << ",serial=" << serial
+       os << class_name() << "(" << "name=" << name << ",serial=" << serial
           << ",hash=" << hashvalue << ",flags=" << flags << ")";
 }
 
@@ -225,11 +224,6 @@ ex symbol::eval(int level) const
 ex symbol::subs(const lst & ls, const lst & lr) const
 {
        GINAC_ASSERT(ls.nops()==lr.nops());
-#ifdef DO_GINAC_ASSERT
-       for (unsigned i=0; i<ls.nops(); i++)
-               GINAC_ASSERT(is_ex_exactly_of_type(ls.op(i),symbol)||
-                            is_ex_of_type(ls.op(i),idx));
-#endif // def DO_GINAC_ASSERT
 
        for (unsigned i=0; i<ls.nops(); i++) {
                if (is_ex_exactly_of_type(ls.op(i),symbol)) {
diff --git a/ginac/tensor.cpp b/ginac/tensor.cpp
new file mode 100644 (file)
index 0000000..28b520f
--- /dev/null
@@ -0,0 +1,375 @@
+/** @file tensor.cpp
+ *
+ *  Implementation of GiNaC's special tensors. */
+
+/*
+ *  GiNaC Copyright (C) 1999-2001 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
+ *  the Free Software Foundation; either version 2 of the License, or
+ *  (at your option) any later version.
+ *
+ *  This program is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with this program; if not, write to the Free Software
+ *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+ */
+
+#include <stdexcept>
+
+#include "tensor.h"
+#include "idx.h"
+#include "indexed.h"
+#include "relational.h"
+#include "numeric.h"
+#include "archive.h"
+#include "utils.h"
+#include "debugmsg.h"
+
+namespace GiNaC {
+
+GINAC_IMPLEMENT_REGISTERED_CLASS(tensor, basic)
+GINAC_IMPLEMENT_REGISTERED_CLASS(tensdelta, tensor)
+GINAC_IMPLEMENT_REGISTERED_CLASS(tensmetric, tensor)
+GINAC_IMPLEMENT_REGISTERED_CLASS(minkmetric, tensmetric)
+GINAC_IMPLEMENT_REGISTERED_CLASS(tensepsilon, tensor)
+
+//////////
+// default constructor, destructor, copy constructor assignment operator and helpers
+//////////
+
+#define DEFAULT_CTORS(classname) \
+classname::classname() : inherited(TINFO_##classname) \
+{ \
+       debugmsg(#classname " default constructor", LOGLEVEL_CONSTRUCT); \
+} \
+void classname::copy(const classname & other) \
+{ \
+       inherited::copy(other); \
+} \
+void classname::destroy(bool call_parent) \
+{ \
+       if (call_parent) \
+               inherited::destroy(call_parent); \
+}
+
+tensor::tensor(unsigned ti) : inherited(ti)
+{
+       debugmsg("tensor constructor from unsigned", LOGLEVEL_CONSTRUCT); \
+}
+
+DEFAULT_CTORS(tensor)
+DEFAULT_CTORS(tensdelta)
+DEFAULT_CTORS(tensmetric)
+DEFAULT_CTORS(tensepsilon)
+
+minkmetric::minkmetric() : pos_sig(false)
+{
+       debugmsg("minkmetric default constructor", LOGLEVEL_CONSTRUCT);
+       tinfo_key = TINFO_minkmetric;
+}
+
+minkmetric::minkmetric(bool ps) : pos_sig(ps)
+{
+       debugmsg("minkmetric constructor from bool", LOGLEVEL_CONSTRUCT);
+       tinfo_key = TINFO_minkmetric;
+}
+
+void minkmetric::copy(const minkmetric & other)
+{
+       inherited::copy(other);
+       pos_sig = other.pos_sig;
+}
+
+void minkmetric::destroy(bool call_parent)
+{
+       if (call_parent)
+               inherited::destroy(call_parent);
+}
+
+//////////
+// archiving
+//////////
+
+#define DEFAULT_ARCHIVING(classname) \
+classname::classname(const archive_node &n, const lst &sym_lst) : inherited(n, sym_lst) \
+{ \
+       debugmsg(#classname " constructor from archive_node", LOGLEVEL_CONSTRUCT); \
+} \
+ex classname::unarchive(const archive_node &n, const lst &sym_lst) \
+{ \
+       return (new classname(n, sym_lst))->setflag(status_flags::dynallocated); \
+} \
+void classname::archive(archive_node &n) const \
+{ \
+       inherited::archive(n); \
+}
+
+DEFAULT_ARCHIVING(tensor)
+DEFAULT_ARCHIVING(tensdelta)
+DEFAULT_ARCHIVING(tensmetric)
+DEFAULT_ARCHIVING(tensepsilon)
+
+minkmetric::minkmetric(const archive_node &n, const lst &sym_lst) : inherited(n, sym_lst)
+{
+       debugmsg("minkmetric constructor from archive_node", LOGLEVEL_CONSTRUCT);
+       n.find_bool("pos_sig", pos_sig);
+}
+
+ex minkmetric::unarchive(const archive_node &n, const lst &sym_lst)
+{
+       return (new minkmetric(n, sym_lst))->setflag(status_flags::dynallocated);
+}
+
+void minkmetric::archive(archive_node &n) const
+{
+       inherited::archive(n);
+       n.add_bool("pos_sig", pos_sig);
+}
+
+//////////
+// functions overriding virtual functions from bases classes
+//////////
+
+#define DEFAULT_COMPARE(classname) \
+int classname::compare_same_type(const basic & other) const \
+{ \
+       /* by default, two tensors of the same class are always identical */ \
+       return 0; \
+}
+
+DEFAULT_COMPARE(tensor)
+DEFAULT_COMPARE(tensdelta)
+DEFAULT_COMPARE(tensmetric)
+DEFAULT_COMPARE(tensepsilon)
+
+int minkmetric::compare_same_type(const basic & other) const
+{
+       GINAC_ASSERT(is_of_type(other, minkmetric));
+       const minkmetric &o = static_cast<const minkmetric &>(other);
+
+       if (pos_sig != o.pos_sig)
+               return pos_sig ? -1 : 1;
+       else
+               return inherited::compare_same_type(other);
+}
+
+void tensdelta::print(std::ostream & os, unsigned upper_precedence) const
+{
+       debugmsg("tensdelta print",LOGLEVEL_PRINT);
+       os << "delta";
+}
+
+void tensmetric::print(std::ostream & os, unsigned upper_precedence) const
+{
+       debugmsg("tensmetric print",LOGLEVEL_PRINT);
+       os << "g";
+}
+
+void minkmetric::print(std::ostream & os, unsigned upper_precedence) const
+{
+       debugmsg("minkmetric print",LOGLEVEL_PRINT);
+       os << "eta";
+}
+
+void tensepsilon::print(std::ostream & os, unsigned upper_precedence) const
+{
+       debugmsg("tensepsilon print",LOGLEVEL_PRINT);
+       os << "eps";
+}
+
+/** Automatic symbolic evaluation of an indexed delta tensor. */
+ex tensdelta::eval_indexed(const basic & i) const
+{
+       GINAC_ASSERT(is_of_type(i, indexed));
+       GINAC_ASSERT(i.nops() == 3);
+       GINAC_ASSERT(is_ex_of_type(i.op(0), tensdelta));
+
+       const idx & i1 = ex_to_idx(i.op(1));
+       const idx & i2 = ex_to_idx(i.op(2));
+
+       // Trace of delta tensor is the dimension of the space
+       if (is_dummy_pair(i1, i2))
+               return i1.get_dim();
+
+       // No further simplifications
+       return i.hold();
+}
+
+/** Automatic symbolic evaluation of an indexed metric tensor. */
+ex tensmetric::eval_indexed(const basic & i) const
+{
+       GINAC_ASSERT(is_of_type(i, indexed));
+       GINAC_ASSERT(i.nops() == 3);
+       GINAC_ASSERT(is_ex_of_type(i.op(0), tensmetric));
+       GINAC_ASSERT(is_ex_of_type(i.op(1), varidx));
+       GINAC_ASSERT(is_ex_of_type(i.op(2), varidx));
+
+       const varidx & i1 = ex_to_varidx(i.op(1));
+       const varidx & i2 = ex_to_varidx(i.op(2));
+
+       // A metric tensor with one covariant and one contravariant index gets
+       // replaced by a delta tensor
+       if (i1.is_covariant() != i2.is_covariant())
+               return delta_tensor(i1, i2);
+
+       // No further simplifications
+       return i.hold();
+}
+
+/** Automatic symbolic evaluation of an indexed Lorentz metric tensor. */
+ex minkmetric::eval_indexed(const basic & i) const
+{
+       GINAC_ASSERT(is_of_type(i, indexed));
+       GINAC_ASSERT(i.nops() == 3);
+       GINAC_ASSERT(is_ex_of_type(i.op(0), minkmetric));
+       GINAC_ASSERT(is_ex_of_type(i.op(1), varidx));
+       GINAC_ASSERT(is_ex_of_type(i.op(2), varidx));
+
+       const varidx & i1 = ex_to_varidx(i.op(1));
+       const varidx & i2 = ex_to_varidx(i.op(2));
+
+       // Numeric evaluation
+       if (static_cast<const indexed &>(i).all_index_values_are(info_flags::nonnegint)) {
+               int n1 = ex_to_numeric(i1.get_value()).to_int(), n2 = ex_to_numeric(i2.get_value()).to_int();
+               if (n1 != n2)
+                       return _ex0();
+               else if (n1 == 0)
+                       return pos_sig ? _ex_1() : _ex1();
+               else
+                       return pos_sig ? _ex1() : _ex_1();
+       }
+
+       // Perform the usual evaluations of a metric tensor
+       return inherited::eval_indexed(i);
+}
+
+/** Contraction of an indexed delta tensor with something else. */
+bool tensdelta::contract_with(ex & self, ex & other) const
+{
+       GINAC_ASSERT(is_ex_of_type(self, indexed));
+       GINAC_ASSERT(is_ex_of_type(other, indexed));
+       GINAC_ASSERT(self.nops() == 3);
+       GINAC_ASSERT(is_ex_of_type(self.op(0), tensdelta));
+
+       // Try to contract first index
+       const idx *self_idx = &ex_to_idx(self.op(1));
+       const idx *free_idx = &ex_to_idx(self.op(2));
+       bool first_index_tried = false;
+
+again:
+       if (self_idx->is_symbolic()) {
+               for (int i=1; i<other.nops(); i++) {
+                       const idx &other_idx = ex_to_idx(other.op(i));
+                       if (is_dummy_pair(*self_idx, other_idx)) {
+
+                               // Contraction found, remove delta tensor and substitute
+                               // index in second object
+                               self = _ex1();
+                               other = other.subs(other_idx == *free_idx);
+                               return true;
+                       }
+               }
+       }
+
+       if (!first_index_tried) {
+
+               // No contraction with first index found, try second index
+               self_idx = &ex_to_idx(self.op(2));
+               free_idx = &ex_to_idx(self.op(1));
+               first_index_tried = true;
+               goto again;
+       }
+
+       return false;
+}
+
+/** Contraction of an indexed metric tensor with something else. */
+bool tensmetric::contract_with(ex & self, ex & other) const
+{
+       GINAC_ASSERT(is_ex_of_type(self, indexed));
+       GINAC_ASSERT(is_ex_of_type(other, indexed));
+       GINAC_ASSERT(self.nops() == 3);
+       GINAC_ASSERT(is_ex_of_type(self.op(0), tensmetric));
+
+       // If contracting with the delta tensor, let the delta do it
+       // (don't raise/lower delta indices)
+       if (is_ex_exactly_of_type(other.op(0), tensdelta))
+               return false;
+
+       // Try to contract first index
+       const idx *self_idx = &ex_to_idx(self.op(1));
+       const idx *free_idx = &ex_to_idx(self.op(2));
+       bool first_index_tried = false;
+
+again:
+       if (self_idx->is_symbolic()) {
+               for (int i=1; i<other.nops(); i++) {
+                       const idx &other_idx = ex_to_idx(other.op(i));
+                       if (is_dummy_pair(*self_idx, other_idx)) {
+
+                               // Contraction found, remove metric tensor and substitute
+                               // index in second object
+                               self = _ex1();
+                               other = other.subs(other_idx == *free_idx);
+                               return true;
+                       }
+               }
+       }
+
+       if (!first_index_tried) {
+
+               // No contraction with first index found, try second index
+               self_idx = &ex_to_idx(self.op(2));
+               free_idx = &ex_to_idx(self.op(1));
+               first_index_tried = true;
+               goto again;
+       }
+
+       return false;
+}
+
+//////////
+// global functions
+//////////
+
+ex delta_tensor(const ex & i1, const ex & i2)
+{
+       if (!is_ex_of_type(i1, idx) || !is_ex_of_type(i2, idx))
+               throw(std::invalid_argument("indices of delta tensor must be of type idx"));
+
+       return indexed(tensdelta(), indexed::symmetric, i1, i2);
+}
+
+ex metric_tensor(const ex & i1, const ex & i2)
+{
+       if (!is_ex_of_type(i1, varidx) || !is_ex_of_type(i2, varidx))
+               throw(std::invalid_argument("indices of metric tensor must be of type varidx"));
+
+       return indexed(tensmetric(), i1, i2);
+}
+
+ex lorentz_g(const ex & i1, const ex & i2, bool pos_sig)
+{
+       if (!is_ex_of_type(i1, varidx) || !is_ex_of_type(i2, varidx))
+               throw(std::invalid_argument("indices of metric tensor must be of type varidx"));
+
+       return indexed(minkmetric(pos_sig), indexed::symmetric, i1, i2);
+}
+
+ex epsilon_tensor(const ex & i1, const ex & i2)
+{
+       if (!is_ex_of_type(i1, idx) || !is_ex_of_type(i2, idx))
+               throw(std::invalid_argument("indices of epsilon tensor must be of type idx"));
+       if (!ex_to_idx(i1).get_dim().is_equal(_ex2()) || !ex_to_idx(i2).get_dim().is_equal(_ex2()))
+               throw(std::invalid_argument("index dimension of epsilon tensor must match number of indices"));
+
+       return indexed(tensepsilon(), indexed::antisymmetric, i1, i2);
+}
+
+} // namespace GiNaC
diff --git a/ginac/tensor.h b/ginac/tensor.h
new file mode 100644 (file)
index 0000000..220b69b
--- /dev/null
@@ -0,0 +1,158 @@
+/** @file tensor.h
+ *
+ *  Interface to GiNaC's special tensors. */
+
+/*
+ *  GiNaC Copyright (C) 1999-2001 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
+ *  the Free Software Foundation; either version 2 of the License, or
+ *  (at your option) any later version.
+ *
+ *  This program is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with this program; if not, write to the Free Software
+ *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+ */
+
+#ifndef __GINAC_TENSOR_H__
+#define __GINAC_TENSOR_H__
+
+#include "ex.h"
+
+namespace GiNaC {
+
+
+/** This class holds one of GiNaC's predefined special tensors such as the
+ *  delta and the metric tensors. They are represented without indices.
+ *  To attach indices to them, wrap them in an object of class indexed. */
+class tensor : public basic
+{
+       GINAC_DECLARE_REGISTERED_CLASS(tensor, basic)
+
+       // other constructors
+protected:
+       tensor(unsigned ti);
+
+       // functions overriding virtual functions from bases classes
+protected:
+       unsigned return_type(void) const { return return_types::noncommutative_composite; }
+};
+
+
+/** This class represents the delta tensor. If indexed, it must have exactly
+ *  two indices of the same type. */
+class tensdelta : public tensor
+{
+       GINAC_DECLARE_REGISTERED_CLASS(tensdelta, tensor)
+
+       // functions overriding virtual functions from bases classes
+public:
+       void print(std::ostream & os, unsigned upper_precedence=0) const;
+       ex eval_indexed(const basic & i) const;
+       bool contract_with(ex & self, ex & other) const;
+};
+
+
+/** This class represents a general metric tensor which can be used to
+ *  raise/lower indices. If indexed, it must have exactly two indices of the
+ *  same type which must be of class varidx or a subclass. */
+class tensmetric : public tensor
+{
+       GINAC_DECLARE_REGISTERED_CLASS(tensmetric, tensor)
+
+       // functions overriding virtual functions from bases classes
+public:
+       void print(std::ostream & os, unsigned upper_precedence=0) const;
+       ex eval_indexed(const basic & i) const;
+       bool contract_with(ex & self, ex & other) const;
+};
+
+
+/** This class represents a Minkowski metric tensor. It has all the
+ *  properties of a metric tensor and is (as a matrix) equal to
+ *  diag(1,-1,-1,...) or diag(-1,1,1,...). */
+class minkmetric : public tensmetric
+{
+       GINAC_DECLARE_REGISTERED_CLASS(minkmetric, tensmetric)
+
+       // other constructors
+public:
+       /** Construct Lorentz metric tensor with given signature. */
+       minkmetric(bool pos_sig);
+
+       // functions overriding virtual functions from bases classes
+public:
+       void print(std::ostream & os, unsigned upper_precedence=0) const;
+       ex eval_indexed(const basic & i) const;
+
+       // member variables
+private:
+       bool pos_sig; /**< If true, the metric is diag(-1,1,1...). Otherwise it is diag(1,-1,-1,...). */
+};
+
+
+/** This class represents the totally antisymmetric epsilon tensor. If
+ *  indexed, all indices must be of the same type and their number must
+ *  be equal to the dimension of the index space. */
+class tensepsilon : public tensor
+{
+       GINAC_DECLARE_REGISTERED_CLASS(tensepsilon, tensor)
+
+       // functions overriding virtual functions from bases classes
+public:
+       void print(std::ostream & os, unsigned upper_precedence=0) const;
+};
+
+
+// utility functions
+inline const tensor &ex_to_tensor(const ex &e)
+{
+       return static_cast<const tensor &>(*e.bp);
+}
+
+/** Create a delta tensor with specified indices. The indices must be of class
+ *  idx or a subclass. The delta tensor is always symmetric and its trace is
+ *  the dimension of the index space.
+ *
+ *  @param i1 First index
+ *  @param i2 Second index
+ *  @return newly constructed delta tensor */
+ex delta_tensor(const ex & i1, const ex & i2);
+
+/** Create a metric tensor with specified indices. The indices must be of
+ *  class varidx or a subclass. A metric tensor with one covariant and one
+ *  contravariant index is equivalent to the delta tensor.
+ *
+ *  @param i1 First index
+ *  @param i2 Second index
+ *  @return newly constructed metric tensor */
+ex metric_tensor(const ex & i1, const ex & i2);
+
+/** Create a Minkowski metric tensor with specified indices. The indices
+ *  must be of class varidx or a subclass. The Lorentz metric is a symmetric
+ *  tensor with a matrix representation of diag(1,-1,-1,...) (negative
+ *  signature, the default) or diag(-1,1,1,...) (positive signature).
+ *
+ *  @param i1 First index
+ *  @param i2 Second index
+ *  @param pos_sig Whether the signature is positive
+ *  @return newly constructed Lorentz metric tensor */
+ex lorentz_g(const ex & i1, const ex & i2, bool pos_sig = false);
+
+/** Create an epsilon tensor with two indices. The indices must be of class
+ *  idx or a subclass, and have a dimension of 2.
+ *
+ *  @param i1 First index
+ *  @param i2 Second index
+ *  @return newly constructed epsilon tensor */
+ex epsilon_tensor(const ex & i1, const ex & i2);
+
+} // namespace GiNaC
+
+#endif // ndef __GINAC_TENSOR_H__
index 927cf79..2b18959 100644 (file)
 
 namespace GiNaC {
 
-const unsigned TINFO_basic        = 0x00000001U;
+const unsigned TINFO_basic         = 0x00000001U;
 
-const unsigned TINFO_expairseq    = 0x00010001U;
-const unsigned TINFO_add          = 0x00011001U;
-const unsigned TINFO_mul          = 0x00011002U;
+const unsigned TINFO_expairseq     = 0x00010001U;
+const unsigned TINFO_add           = 0x00011001U;
+const unsigned TINFO_mul           = 0x00011002U;
 
-const unsigned TINFO_symbol       = 0x00020001U;
-const unsigned TINFO_constant     = 0x00021001U;
+const unsigned TINFO_symbol        = 0x00020001U;
+const unsigned TINFO_constant      = 0x00021001U;
 
-const unsigned TINFO_exprseq      = 0x00030001U;
-const unsigned TINFO_function     = 0x00031001U;
-const unsigned TINFO_ncmul        = 0x00031002U;
+const unsigned TINFO_exprseq       = 0x00030001U;
+const unsigned TINFO_function      = 0x00031001U;
+const unsigned TINFO_ncmul         = 0x00031002U;
 
-const unsigned TINFO_lst          = 0x00040001U;
+const unsigned TINFO_lst           = 0x00040001U;
 
-const unsigned TINFO_matrix       = 0x00050001U;
+const unsigned TINFO_matrix        = 0x00050001U;
 
-const unsigned TINFO_power        = 0x00060001U;
+const unsigned TINFO_power         = 0x00060001U;
 
-const unsigned TINFO_relational   = 0x00070001U;
+const unsigned TINFO_relational    = 0x00070001U;
 
-const unsigned TINFO_fail         = 0x00080001U;
+const unsigned TINFO_fail          = 0x00080001U;
 
-const unsigned TINFO_numeric      = 0x00090001U;
+const unsigned TINFO_numeric       = 0x00090001U;
 
-const unsigned TINFO_pseries      = 0x000a0001U;
+const unsigned TINFO_pseries       = 0x000a0001U;
 
-const unsigned TINFO_indexed      = 0x000b0001U;
-const unsigned TINFO_algebra      = 0x000b1001U;
-const unsigned TINFO_clifford     = 0x000b1101U;
-const unsigned TINFO_color        = 0x000b1201U;
-const unsigned TINFO_isospin      = 0x000b1301U;
-const unsigned TINFO_simp_lor     = 0x000b1401U;
-const unsigned TINFO_lortensor    = 0x000b1501U;
+const unsigned TINFO_indexed       = 0x000b0001U;
 
-const unsigned TINFO_structure    = 0x000c0001U;
-// reserved up to                   0x000cffffU
+const unsigned TINFO_structure     = 0x000c0001U;
+// reserved up to                    0x000cffffU
 // for user defined structures
 
-const unsigned TINFO_idx          = 0x000d0001U;
-const unsigned TINFO_coloridx     = 0x000d1001U;
-const unsigned TINFO_lorentzidx   = 0x000d1002U;
+const unsigned TINFO_idx           = 0x000d0001U;
+const unsigned TINFO_varidx        = 0x000d1001U;
+
+const unsigned TINFO_tensor        = 0x000e0001U;
+const unsigned TINFO_tensdelta     = 0x000e1001U;
+const unsigned TINFO_tensmetric    = 0x000e1002U;
+const unsigned TINFO_minkmetric    = 0x000e2001U;
+const unsigned TINFO_tensepsilon   = 0x000e1003U;
 
 } // namespace GiNaC