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
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) {
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;
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
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
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
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
}
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;
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
{
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 << ",";
#include "numeric.h"
#include "power.h"
#include "relational.h"
+#include "indexed.h"
#include "input_lexer.h"
#include "debugmsg.h"
#include "utils.h"
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
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
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;
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);
class expand_options {
public:
enum {
- expand_trigonometric = 0x0001
+ expand_trigonometric = 0x0001,
+ expand_indexed = 0x0002
};
};
has_indices, // object has at least one index
// answered by class idx
- idx,
-
- // answered by class coloridx
- coloridx,
-
- // answered by class lorentzidx
- lorentzidx
+ idx
};
};
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);
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 << ")"
#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__
#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 << ")";
}
{
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
#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:
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
/** @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 {
// 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());
}
//////////
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. */
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);
}
//////////
// 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
/** @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__
#include "archive.h"
#include "numeric.h"
#include "lst.h"
+#include "idx.h"
+#include "indexed.h"
#include "utils.h"
#include "debugmsg.h"
#include "power.h"
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)
return false;
}
-/** evaluate matrix entry by entry. */
+/** Evaluate matrix entry by entry. */
ex matrix::eval(int level) const
{
debugmsg("matrix eval",LOGLEVEL_MEMBER_FUNCTION);
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);
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
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
//////////
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
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!"));
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;
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));
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:
// 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) << ")";
}
{
debugmsg("power printraw",LOGLEVEL_PRINT);
- os << "power(";
+ os << class_name() << "(";
basis.printraw(os);
os << ",";
exponent.printraw(os);
{
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);
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;
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 << ")";
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;
++it;
++o_it;
}
+
// so they are equal.
return 0;
}
void relational::printraw(std::ostream & os) const
{
debugmsg("relational printraw",LOGLEVEL_PRINT);
- os << "RELATIONAL(";
+ os << class_name() << "(";
lh.printraw(os);
os << ",";
rh.printraw(os);
{
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
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
#include "symbol.h"
#include "lst.h"
-#include "idx.h"
#include "archive.h"
#include "debugmsg.h"
#include "utils.h"
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 << ")";
}
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)) {
--- /dev/null
+/** @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
--- /dev/null
+/** @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__
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
git://www.ginac.de / ginac.git / commitdiff