From 094911eb78cacb6f2877a70c9ac74766df58ccea Mon Sep 17 00:00:00 2001 From: Christian Bauer Date: Tue, 6 Mar 2001 20:07:16 +0000 Subject: [PATCH] - revamped indexed objects - subs() works on matrices --- ginac/Makefile.am | 19 +- ginac/add.cpp | 11 - ginac/add.h | 2 +- ginac/basic.cpp | 32 +- ginac/basic.h | 4 +- ginac/constant.cpp | 2 +- ginac/container.pl | 2 +- ginac/ex.cpp | 28 +- ginac/ex.h | 5 +- ginac/expairseq.cpp | 3 +- ginac/flags.h | 11 +- ginac/function.pl | 4 +- ginac/ginac.h | 10 +- ginac/idx.cpp | 540 ++++++++++--------------------- ginac/idx.h | 116 ++++--- ginac/indexed.cpp | 733 +++++++++++++++++++++++++++++++++++-------- ginac/indexed.h | 204 ++++++++++-- ginac/matrix.cpp | 201 +++++++++++- ginac/matrix.h | 4 +- ginac/mul.cpp | 13 - ginac/mul.h | 2 +- ginac/ncmul.cpp | 14 - ginac/ncmul.h | 2 +- ginac/numeric.cpp | 2 +- ginac/power.cpp | 6 +- ginac/power.h | 1 + ginac/pseries.cpp | 5 +- ginac/relational.cpp | 2 +- ginac/structure.cpp | 4 +- ginac/structure.pl | 4 +- ginac/symbol.cpp | 8 +- ginac/tensor.cpp | 375 ++++++++++++++++++++++ ginac/tensor.h | 158 ++++++++++ ginac/tinfos.h | 55 ++-- 34 files changed, 1887 insertions(+), 695 deletions(-) create mode 100644 ginac/tensor.cpp create mode 100644 ginac/tensor.h diff --git a/ginac/Makefile.am b/ginac/Makefile.am index 7a3cd72d..fffe951c 100644 --- a/ginac/Makefile.am +++ b/ginac/Makefile.am @@ -5,20 +5,17 @@ libginac_la_SOURCES = add.cpp archive.cpp basic.cpp constant.cpp ex.cpp \ expairseq.cpp exprseq.cpp fail.cpp function.cpp inifcns.cpp \ inifcns_trans.cpp inifcns_zeta.cpp inifcns_gamma.cpp matrix.cpp mul.cpp \ normal.cpp numeric.cpp operators.cpp power.cpp registrar.cpp relational.cpp \ - symbol.cpp pseries.cpp utils.cpp ncmul.cpp clifford.cpp structure.cpp \ - color.cpp indexed.cpp idx.cpp isospin.cpp exprseq_suppl.cpp lst.cpp \ - lst_suppl.cpp simp_lor.cpp coloridx.cpp lorentzidx.cpp lortensor.cpp \ - input_parser.yy input_lexer.ll input_lexer.h \ - remember.h remember.cpp debugmsg.h utils.h + symbol.cpp pseries.cpp utils.cpp ncmul.cpp structure.cpp exprseq_suppl.cpp \ + lst.cpp lst_suppl.cpp input_parser.yy input_lexer.ll input_lexer.h \ + remember.h remember.cpp debugmsg.h utils.h idx.cpp indexed.cpp tensor.cpp libginac_la_LDFLAGS = -version-info $(LT_CURRENT):$(LT_REVISION):$(LT_AGE) \ -release $(LT_RELEASE) ginacincludedir = $(includedir)/ginac -ginacinclude_HEADERS = ginac.h add.h archive.h basic.h clifford.h color.h \ - coloridx.h constant.h ex.h expair.h expairseq.h exprseq.h fail.h flags.h \ - function.h idx.h indexed.h inifcns.h isospin.h lorentzidx.h lst.h matrix.h \ - mul.h ncmul.h normal.h numeric.h operators.h power.h registrar.h \ - relational.h pseries.h simp_lor.h structure.h symbol.h lortensor.h \ - tinfos.h assertion.h version.h +ginacinclude_HEADERS = ginac.h add.h archive.h basic.h constant.h ex.h \ + expair.h expairseq.h exprseq.h fail.h flags.h function.h inifcns.h \ + lst.h matrix.h mul.h ncmul.h normal.h numeric.h operators.h power.h \ + registrar.h relational.h pseries.h structure.h symbol.h tinfos.h assertion.h \ + version.h idx.h indexed.h tensor.h LFLAGS = -Pginac_yy -olex.yy.c YFLAGS = -p ginac_yy -d EXTRA_DIST = container.pl function.pl structure.pl input_parser.h version.h.in diff --git a/ginac/add.cpp b/ginac/add.cpp index 0c22465a..6c6a1638 100644 --- a/ginac/add.cpp +++ b/ginac/add.cpp @@ -354,17 +354,6 @@ ex add::eval(int level) const return this->hold(); } -exvector add::get_indices(void) const -{ - // FIXME: all terms in the sum should have the same indices (compatible - // tensors) however this is not checked, since there is no function yet - // which compares indices (idxvector can be unsorted) - if (seq.size()==0) { - return exvector(); - } - return (seq.begin())->rest.get_indices(); -} - ex add::simplify_ncmul(const exvector & v) const { if (seq.size()==0) { diff --git a/ginac/add.h b/ginac/add.h index 8a93186b..8738bdf3 100644 --- a/ginac/add.h +++ b/ginac/add.h @@ -59,7 +59,7 @@ public: numeric integer_content(void) const; ex smod(const numeric &xi) const; numeric max_coefficient(void) const; - exvector get_indices(void) const; + exvector get_free_indices(void) const; ex simplify_ncmul(const exvector & v) const; protected: ex derivative(const symbol & s) const; diff --git a/ginac/basic.cpp b/ginac/basic.cpp index 6adc4c49..2414c929 100644 --- a/ginac/basic.cpp +++ b/ginac/basic.cpp @@ -117,7 +117,7 @@ void basic::archive(archive_node &n) const void basic::print(std::ostream & os, unsigned upper_precedence) const { debugmsg("basic print",LOGLEVEL_PRINT); - os << "[basic object]"; + os << "[" << class_name() << " object]"; } /** Output to ostream in ugly raw format, so brave developers can have a look @@ -125,7 +125,7 @@ void basic::print(std::ostream & os, unsigned upper_precedence) const void basic::printraw(std::ostream & os) const { debugmsg("basic printraw",LOGLEVEL_PRINT); - os << "[basic object]"; + os << "[" << class_name() << " object]"; } /** Output to ostream formatted in tree- (indented-) form, so developers can @@ -287,6 +287,31 @@ ex basic::evalf(int level) const return *this; } +/** Perform automatic symbolic evaluations on indexed expression that + * contains this object as the base expression. */ +ex basic::eval_indexed(const basic & i) const + // this function can't take a "const ex & i" because that would result + // in an infinite eval() loop +{ + // There is nothing to do for basic objects + return i.hold(); +} + +/** Try to contract two indexed expressions. If a contraction exists, the + * function overwrites one or both arguments and returns true. Otherwise it + * returns false. It is guaranteed that both expressions are of class + * indexed (or a subclass) and that at least one dummy index has been + * found. + * + * @param self The first indexed expression; it's base object is *this + * @param other The second indexed expression + * @return true if the contraction was successful, false otherwise */ +bool basic::contract_with(ex & self, ex & other) const +{ + // Do nothing + return false; +} + /** Substitute a set of symbols by arbitrary expressions. The ex returned * will already be evaluated. */ ex basic::subs(const lst & ls, const lst & lr) const @@ -320,7 +345,8 @@ ex basic::diff(const symbol & s, unsigned nth) const return ndiff; } -exvector basic::get_indices(void) const +/** Return a vector containing the free indices of an expression. */ +exvector basic::get_free_indices(void) const { return exvector(); // return an empty exvector } diff --git a/ginac/basic.h b/ginac/basic.h index e1325745..842417ca 100644 --- a/ginac/basic.h +++ b/ginac/basic.h @@ -127,8 +127,10 @@ public: // only const functions please (may break reference counting) virtual numeric integer_content(void) const; virtual ex smod(const numeric &xi) const; virtual numeric max_coefficient(void) const; - virtual exvector get_indices(void) const; + virtual exvector get_free_indices(void) const; virtual ex simplify_ncmul(const exvector & v) const; + virtual ex eval_indexed(const basic & i) const; + virtual bool contract_with(ex & self, ex & other) const; protected: // non-const functions should be called from class ex only virtual ex derivative(const symbol & s) const; virtual int compare_same_type(const basic & other) const; diff --git a/ginac/constant.cpp b/ginac/constant.cpp index 56f5a851..90025154 100644 --- a/ginac/constant.cpp +++ b/ginac/constant.cpp @@ -138,7 +138,7 @@ void constant::print(std::ostream & os, unsigned upper_precedence) const void constant::printraw(std::ostream & os) const { debugmsg("constant printraw",LOGLEVEL_PRINT); - os << "constant(" << name << ")"; + os << class_name() << "(" << name << ")"; } void constant::printtree(std::ostream & os, unsigned indent) const diff --git a/ginac/container.pl b/ginac/container.pl index 80fb9fac..571330f9 100755 --- a/ginac/container.pl +++ b/ginac/container.pl @@ -410,7 +410,7 @@ void ${CONTAINER}::printraw(std::ostream & os) const { debugmsg("${CONTAINER} printraw",LOGLEVEL_PRINT); - os << "${CONTAINER}("; + os << class_name() << "("; for (${STLT}::const_iterator cit=seq.begin(); cit!=seq.end(); ++cit) { (*cit).bp->printraw(os); os << ","; diff --git a/ginac/ex.cpp b/ginac/ex.cpp index 6acde9f9..55482551 100644 --- a/ginac/ex.cpp +++ b/ginac/ex.cpp @@ -37,6 +37,7 @@ #include "numeric.h" #include "power.h" #include "relational.h" +#include "indexed.h" #include "input_lexer.h" #include "debugmsg.h" #include "utils.h" @@ -241,10 +242,33 @@ ex ex::subs(const ex & e) const return bp->subs(e); } -exvector ex::get_indices(void) const +/** Return a vector containing the free indices of the object. */ +exvector ex::get_free_indices(void) const { GINAC_ASSERT(bp!=0); - return bp->get_indices(); + return bp->get_free_indices(); +} + +/** Simplify/canonicalize expression containing indexed objects. This + * performs contraction of dummy indices where possible and checks whether + * the free indices in sums are consistent. + * + * @return simplified expression */ +ex ex::simplify_indexed(void) const +{ + return GiNaC::simplify_indexed(*this); +} + +/** Simplify/canonicalize expression containing indexed objects. This + * performs contraction of dummy indices where possible, checks whether + * the free indices in sums are consistent, and automatically replaces + * scalar products by known values if desired. + * + * @param sp Scalar products to be replaced automatically + * @return simplified expression */ +ex ex::simplify_indexed(const scalar_products & sp) const +{ + return GiNaC::simplify_indexed(*this, sp); } ex ex::simplify_ncmul(const exvector & v) const diff --git a/ginac/ex.h b/ginac/ex.h index 0e2adaf0..f2feba73 100644 --- a/ginac/ex.h +++ b/ginac/ex.h @@ -35,6 +35,7 @@ extern const ex & _ex0(void); ///< single ex(numeric(0)) class symbol; class lst; +class scalar_products; /** Lightweight wrapper for GiNaC's symbolic objects. Basically all it does is * to hold a pointer to the other objects, manage the reference counting and @@ -107,7 +108,9 @@ public: ex series(const ex & r, int order, unsigned options = 0) const; ex subs(const lst & ls, const lst & lr) const; ex subs(const ex & e) const; - exvector get_indices(void) const; + exvector get_free_indices(void) const; + ex simplify_indexed(void) const; + ex simplify_indexed(const scalar_products & sp) const; ex simplify_ncmul(const exvector & v) const; ex operator[](const ex & index) const; ex operator[](int i) const; diff --git a/ginac/expairseq.cpp b/ginac/expairseq.cpp index bad422cd..cfcfa16d 100644 --- a/ginac/expairseq.cpp +++ b/ginac/expairseq.cpp @@ -207,8 +207,7 @@ void expairseq::print(std::ostream &os, unsigned upper_precedence) const void expairseq::printraw(std::ostream &os) const { debugmsg("expairseq printraw",LOGLEVEL_PRINT); - - os << "expairseq("; + os << class_name() << "("; for (epvector::const_iterator cit=seq.begin(); cit!=seq.end(); ++cit) { os << "("; (*cit).rest.printraw(os); diff --git a/ginac/flags.h b/ginac/flags.h index 93506f10..987a4df0 100644 --- a/ginac/flags.h +++ b/ginac/flags.h @@ -28,7 +28,8 @@ namespace GiNaC { class expand_options { public: enum { - expand_trigonometric = 0x0001 + expand_trigonometric = 0x0001, + expand_indexed = 0x0002 }; }; @@ -128,13 +129,7 @@ public: has_indices, // object has at least one index // answered by class idx - idx, - - // answered by class coloridx - coloridx, - - // answered by class lorentzidx - lorentzidx + idx }; }; diff --git a/ginac/function.pl b/ginac/function.pl index 0a157661..ae9c832f 100755 --- a/ginac/function.pl +++ b/ginac/function.pl @@ -642,7 +642,7 @@ void function::printraw(std::ostream & os) const GINAC_ASSERT(serialprint(os); @@ -666,7 +666,7 @@ void function::printtree(std::ostream & os, unsigned indent) const GINAC_ASSERT(serial #include "idx.h" -#include "ex.h" +#include "symbol.h" #include "lst.h" -#include "relational.h" #include "archive.h" #include "utils.h" #include "debugmsg.h" @@ -33,160 +32,118 @@ 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; isetflag(status_flags::dynallocated); } -/** Archive the object. */ void idx::archive(archive_node &n) const { inherited::archive(n); - n.add_bool("symbolic", symbolic); + n.add_ex("value", value); + n.add_ex("dim", dim); +} + +void varidx::archive(archive_node &n) const +{ + inherited::archive(n); n.add_bool("covariant", covariant); - if (symbolic) - n.add_string("name", name); - else - n.add_unsigned("value", value); } ////////// // functions overriding virtual functions from bases classes ////////// -// public - void idx::printraw(std::ostream & os) const { - debugmsg("idx printraw",LOGLEVEL_PRINT); - - os << "idx("; - - if (symbolic) { - os << "symbolic,name=" << name; - } else { - os << "non symbolic,value=" << value; - } + debugmsg("idx printraw", LOGLEVEL_PRINT); - if (covariant) { - os << ",covariant"; - } else { - os << ",contravariant"; - } - - os << ",serial=" << serial; + os << class_name() << "("; + value.printraw(os); + os << ",dim="; + dim.printraw(os); os << ",hash=" << hashvalue << ",flags=" << flags; os << ")"; } @@ -195,334 +152,173 @@ void idx::printtree(std::ostream & os, unsigned indent) const { debugmsg("idx printtree",LOGLEVEL_PRINT); - os << std::string(indent,' ') << "idx: "; - - if (symbolic) { - os << "symbolic,name=" << name; - } else { - os << "non symbolic,value=" << value; - } - - if (covariant) { - os << ",covariant"; - } else { - os << ",contravariant"; - } - - os << ", serial=" << serial - << ", hash=" << hashvalue + os << std::string(indent, ' ') << "type=" << class_name(); + value.printtree(os, indent + delta_indent); + os << std::string(indent, ' '); + os << ", hash=" << hashvalue << " (0x" << std::hex << hashvalue << std::dec << ")" << ", flags=" << flags << std::endl; } void idx::print(std::ostream & os, unsigned upper_precedence) const { - debugmsg("idx print",LOGLEVEL_PRINT); + debugmsg("idx print", LOGLEVEL_PRINT); - if (covariant) { - os << "_"; - } else { - os << "~"; - } - if (symbolic) { - os << name; - } else { - os << value; - } -} + os << "."; -bool idx::info(unsigned inf) const -{ - if (inf==info_flags::idx) return true; - return inherited::info(inf); + bool need_parens = !(is_ex_exactly_of_type(value, numeric) || is_ex_of_type(value, symbol)); + if (need_parens) + os << "("; + os << value; + if (need_parens) + os << ")"; } -ex idx::subs(const lst & ls, const lst & lr) const +void varidx::print(std::ostream & os, unsigned upper_precedence) const { - GINAC_ASSERT(ls.nops()==lr.nops()); -#ifdef DO_GINAC_ASSERT - for (unsigned i=0; i(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(other); - if (symbolic && o.symbolic) { - // both symbolic: compare serials - if (serial==o.serial) { - return 0; - } - return serial(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; iindex + 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(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_cast(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(duplicate()); - i_copy->covariant = !i_copy->covariant; - i_copy->clearflag(status_flags::hash_calculated); + idx *i_copy = static_cast(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(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(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(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 -//#include -#include "basic.h" #include "ex.h" namespace GiNaC { -/** This class holds one index of an indexed object. Indices can be symbolic - * (e.g. "mu", "i") or numeric (unsigned integer), and they can be contravariant - * (the default) or covariant. */ +/** This class holds one index of an indexed object. Indices can + * theoretically consist of any symbolic expression but they are usually + * only just a symbol (e.g. "mu", "i") or numeric (integer). Indices belong + * to a space with a certain numeric or symbolic dimension. */ class idx : public basic { GINAC_DECLARE_REGISTERED_CLASS(idx, basic) // other constructors public: - explicit idx(bool cov); - explicit idx(const std::string & n, bool cov=false); - explicit idx(const char * n, bool cov=false); - explicit idx(unsigned v, bool cov=false); + /** Construct index with given value and dimension. + * + * @param v Value of index (numeric or symbolic) + * @param dim Dimension of index space (numeric or symbolic) + * @return newly constructed index */ + explicit idx(const ex & v, const ex & dim); // functions overriding virtual functions from bases classes public: @@ -52,56 +52,96 @@ public: void print(std::ostream & os, unsigned upper_precedence=0) const; bool info(unsigned inf) const; protected: - bool is_equal_same_type(const basic & other) const; - unsigned calchash(void) const; ex subs(const lst & ls, const lst & lr) const; - // new virtual functions which can be overridden by derived classes + // new virtual functions in this class public: - virtual bool is_co_contra_pair(const basic & other) const; - virtual ex toggle_covariant(void) const; + /** Check whether the index forms a dummy index pair with another index + * of the same type. */ + virtual bool is_dummy_pair_same_type(const basic & other) const; // non-virtual functions in this class public: - /** Check whether index is symbolic (not numeric). */ - bool is_symbolic(void) const {return symbolic;} + /** Get value of index. */ + ex get_value(void) const {return value;} - /** Get numeric value of index. Undefined for symbolic indices. */ - unsigned get_value(void) const {return value;} + /** Check whether the index is numeric. */ + bool is_numeric(void) const {return is_ex_exactly_of_type(value, numeric);} - /** Check whether index is covariant (not contravariant). */ + /** Check whether the index is symbolic. */ + bool is_symbolic(void) const {return !is_ex_exactly_of_type(value, numeric);} + + /** Get dimension of index space. */ + ex get_dim(void) const {return dim;} + + /** Check whether the dimension is numeric. */ + bool is_dim_numeric(void) const {return is_ex_exactly_of_type(dim, numeric);} + + /** Check whether the dimension is symbolic. */ + bool is_dim_symbolic(void) const {return !is_ex_exactly_of_type(dim, numeric);} + + // member variables +protected: + ex value; /**< Expression that constitutes the index (numeric or symbolic name) */ + ex dim; /**< Dimension of space (can be symbolic or numeric) */ +}; + + +/** This class holds an index with a variance (co- or contravariant). There + * is an associated metric tensor that can be used to raise/lower indices. */ +class varidx : public idx +{ + GINAC_DECLARE_REGISTERED_CLASS(varidx, idx) + + // other constructors +public: + /** Construct index with given value, dimension and variance. + * + * @param v Value of index (numeric or symbolic) + * @param dim Dimension of index space (numeric or symbolic) + * @param covariant Make covariant index (default is contravariant) + * @return newly constructed index */ + varidx(const ex & v, const ex & dim, bool covariant = false); + + // functions overriding virtual functions from bases classes +public: + void print(std::ostream & os, unsigned upper_precedence=0) const; + bool is_dummy_pair_same_type(const basic & other) const; + + // non-virtual functions in this class +public: + /** Check whether the index is covariant. */ bool is_covariant(void) const {return covariant;} - void setname(const std::string & n) {name=n;} - std::string getname(void) const {return name;} + /** Check whether the index is contravariant (not covariant). */ + bool is_contravariant(void) const {return !covariant;} -private: - std::string & autoname_prefix(void); + /** Make a new index with the same value but the opposite variance. */ + ex toggle_variance(void) const; // member variables protected: - unsigned serial; /**< Unique serial number for comparing symbolic indices */ - bool symbolic; /**< Is index symbolic? */ - std::string name; /**< Symbolic name (if symbolic == true) */ - unsigned value; /**< Numeric value (if symbolic == false) */ - static unsigned next_serial; - bool covariant; /**< x_mu, default is contravariant: x~mu */ + bool covariant; /**< x.mu, default is contravariant: x~mu */ }; + // utility functions -inline const idx &ex_to_idx(const ex &e) +inline const idx &ex_to_idx(const ex & e) { return static_cast(*e.bp); } -// global functions +inline const varidx &ex_to_varidx(const ex & e) +{ + return static_cast(*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 diff --git a/ginac/indexed.cpp b/ginac/indexed.cpp index 9456f38c..ea4266cb 100644 --- a/ginac/indexed.cpp +++ b/ginac/indexed.cpp @@ -1,6 +1,6 @@ /** @file indexed.cpp * - * Implementation of GiNaC's index carrying objects. */ + * Implementation of GiNaC's indexed expressions. */ /* * GiNaC Copyright (C) 1999-2001 Johannes Gutenberg University Mainz, Germany @@ -20,11 +20,16 @@ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */ -#include +#include #include "indexed.h" -#include "ex.h" #include "idx.h" +#include "add.h" +#include "mul.h" +#include "ncmul.h" +#include "power.h" +#include "archive.h" +#include "utils.h" #include "debugmsg.h" namespace GiNaC { @@ -35,107 +40,119 @@ GINAC_IMPLEMENT_REGISTERED_CLASS(indexed, exprseq) // default constructor, destructor, copy constructor assignment operator and helpers ////////// -// public - -indexed::indexed() +indexed::indexed() : symmetry(unknown) { - debugmsg("indexed default constructor",LOGLEVEL_CONSTRUCT); - tinfo_key=TINFO_indexed; + debugmsg("indexed default constructor", LOGLEVEL_CONSTRUCT); + tinfo_key = TINFO_indexed; } -// protected - void indexed::copy(const indexed & other) { inherited::copy(other); + symmetry = other.symmetry; } void indexed::destroy(bool call_parent) { - if (call_parent) { + if (call_parent) inherited::destroy(call_parent); - } } ////////// // other constructors ////////// -// public +indexed::indexed(const ex & b) : inherited(b), symmetry(unknown) +{ + debugmsg("indexed constructor from ex", LOGLEVEL_CONSTRUCT); + tinfo_key = TINFO_indexed; + GINAC_ASSERT(all_indices_of_type_idx()); +} -/** Construct indexed object with one index. The index must be of class idx - * or a subclass. - * - * @param i1 The index - * @return newly constructed indexed object */ -indexed::indexed(const ex & i1) : inherited(i1) +indexed::indexed(const ex & b, const ex & i1) : inherited(b, i1), symmetry(unknown) { - debugmsg("indexed constructor from ex",LOGLEVEL_CONSTRUCT); - tinfo_key=TINFO_indexed; - GINAC_ASSERT(all_of_type_idx()); + debugmsg("indexed constructor from ex,ex", LOGLEVEL_CONSTRUCT); + tinfo_key = TINFO_indexed; + GINAC_ASSERT(all_indices_of_type_idx()); } -/** Construct indexed object with two indices. The indices must be of class - * idx or a subclass. - * - * @param i1 First index - * @param i2 Second index - * @return newly constructed indexed object */ -indexed::indexed(const ex & i1, const ex & i2) : inherited(i1,i2) +indexed::indexed(const ex & b, const ex & i1, const ex & i2) : inherited(b, i1, i2), symmetry(unknown) { - debugmsg("indexed constructor from ex,ex",LOGLEVEL_CONSTRUCT); - tinfo_key=TINFO_indexed; - GINAC_ASSERT(all_of_type_idx()); + debugmsg("indexed constructor from ex,ex,ex", LOGLEVEL_CONSTRUCT); + tinfo_key = TINFO_indexed; + GINAC_ASSERT(all_indices_of_type_idx()); } -/** Construct indexed object with three indices. The indices must be of class - * idx or a subclass. - * - * @param i1 First index - * @param i2 Second index - * @param i3 Third index - * @return newly constructed indexed object */ -indexed::indexed(const ex & i1, const ex & i2, const ex & i3) - : inherited(i1,i2,i3) +indexed::indexed(const ex & b, const ex & i1, const ex & i2, const ex & i3) : inherited(b, i1, i2, i3), symmetry(unknown) { - debugmsg("indexed constructor from ex,ex,ex",LOGLEVEL_CONSTRUCT); - tinfo_key=TINFO_indexed; - GINAC_ASSERT(all_of_type_idx()); + debugmsg("indexed constructor from ex,ex,ex,ex", LOGLEVEL_CONSTRUCT); + tinfo_key = TINFO_indexed; + GINAC_ASSERT(all_indices_of_type_idx()); } -/** Construct indexed object with four indices. The indices must be of class - * idx or a subclass. - * - * @param i1 First index - * @param i2 Second index - * @param i3 Third index - * @param i4 Fourth index - * @return newly constructed indexed object */ -indexed::indexed(const ex & i1, const ex & i2, const ex & i3, const ex & i4) - : inherited(i1,i2,i3,i4) +indexed::indexed(const ex & b, const ex & i1, const ex & i2, const ex & i3, const ex & i4) : inherited(b, i1, i2, i3, i4), symmetry(unknown) +{ + debugmsg("indexed constructor from ex,ex,ex,ex,ex", LOGLEVEL_CONSTRUCT); + tinfo_key = TINFO_indexed; + GINAC_ASSERT(all_indices_of_type_idx()); +} + +indexed::indexed(const ex & b, symmetry_type symm, const ex & i1, const ex & i2) : inherited(b, i1, i2), symmetry(symm) { - debugmsg("indexed constructor from ex,ex,ex,ex",LOGLEVEL_CONSTRUCT); - tinfo_key=TINFO_indexed; - GINAC_ASSERT(all_of_type_idx()); + debugmsg("indexed constructor from ex,symmetry,ex,ex", LOGLEVEL_CONSTRUCT); + tinfo_key = TINFO_indexed; + GINAC_ASSERT(all_indices_of_type_idx()); } -/** Construct indexed object with a specified vector of indices. The indices - * must be of class idx or a subclass. - * - * @param iv Vector of indices - * @return newly constructed indexed object */ -indexed::indexed(const exvector & iv) : inherited(iv) +indexed::indexed(const ex & b, symmetry_type symm, const ex & i1, const ex & i2, const ex & i3) : inherited(b, i1, i2, i3), symmetry(symm) { - debugmsg("indexed constructor from exvector",LOGLEVEL_CONSTRUCT); - tinfo_key=TINFO_indexed; - GINAC_ASSERT(all_of_type_idx()); + debugmsg("indexed constructor from ex,symmetry,ex,ex,ex", LOGLEVEL_CONSTRUCT); + tinfo_key = TINFO_indexed; + GINAC_ASSERT(all_indices_of_type_idx()); +} + +indexed::indexed(const ex & b, symmetry_type symm, const ex & i1, const ex & i2, const ex & i3, const ex & i4) : inherited(b, i1, i2, i3, i4), symmetry(symm) +{ + debugmsg("indexed constructor from ex,symmetry,ex,ex,ex,ex", LOGLEVEL_CONSTRUCT); + tinfo_key = TINFO_indexed; + GINAC_ASSERT(all_indices_of_type_idx()); +} + +indexed::indexed(const ex & b, const exvector & v) : inherited(b), symmetry(unknown) +{ + debugmsg("indexed constructor from ex,exvector", LOGLEVEL_CONSTRUCT); + seq.insert(seq.end(), v.begin(), v.end()); + tinfo_key = TINFO_indexed; + GINAC_ASSERT(all_indices_of_type_idx()); +} + +indexed::indexed(const ex & b, symmetry_type symm, const exvector & v) : inherited(b), symmetry(symm) +{ + debugmsg("indexed constructor from ex,symmetry,exvector", LOGLEVEL_CONSTRUCT); + seq.insert(seq.end(), v.begin(), v.end()); + tinfo_key = TINFO_indexed; + GINAC_ASSERT(all_indices_of_type_idx()); +} + +indexed::indexed(symmetry_type symm, const exprseq & es) : inherited(es), symmetry(symm) +{ + debugmsg("indexed constructor from symmetry,exprseq", LOGLEVEL_CONSTRUCT); + tinfo_key = TINFO_indexed; + GINAC_ASSERT(all_indices_of_type_idx()); } -indexed::indexed(exvector * ivp) : inherited(ivp) +indexed::indexed(symmetry_type symm, const exvector & v, bool discardable) : inherited(v, discardable), symmetry(symm) { - debugmsg("indexed constructor from exvector *",LOGLEVEL_CONSTRUCT); - tinfo_key=TINFO_indexed; - GINAC_ASSERT(all_of_type_idx()); + debugmsg("indexed constructor from symmetry,exvector", LOGLEVEL_CONSTRUCT); + tinfo_key = TINFO_indexed; + GINAC_ASSERT(all_indices_of_type_idx()); +} + +indexed::indexed(symmetry_type symm, exvector * vp) : inherited(vp), symmetry(symm) +{ + debugmsg("indexed constructor from symmetry,exvector *", LOGLEVEL_CONSTRUCT); + tinfo_key = TINFO_indexed; + GINAC_ASSERT(all_indices_of_type_idx()); } ////////// @@ -146,7 +163,9 @@ indexed::indexed(exvector * ivp) : inherited(ivp) indexed::indexed(const archive_node &n, const lst &sym_lst) : inherited(n, sym_lst) { debugmsg("indexed constructor from archive_node", LOGLEVEL_CONSTRUCT); - tinfo_key = TINFO_indexed; + unsigned int symm; + if (!(n.find_unsigned("symmetry", symm))) + throw (std::runtime_error("unknown indexed symmetry type in archive")); } /** Unarchive the object. */ @@ -159,90 +178,202 @@ ex indexed::unarchive(const archive_node &n, const lst &sym_lst) void indexed::archive(archive_node &n) const { inherited::archive(n); + n.add_unsigned("symmetry", symmetry); } ////////// // functions overriding virtual functions from bases classes ////////// -// public - void indexed::printraw(std::ostream & os) const { - debugmsg("indexed printraw",LOGLEVEL_PRINT); - os << "indexed(indices="; + debugmsg("indexed printraw", LOGLEVEL_PRINT); + GINAC_ASSERT(seq.size() > 0); + + os << class_name() << "("; + seq[0].printraw(os); + os << ",indices="; printrawindices(os); os << ",hash=" << hashvalue << ",flags=" << flags << ")"; } void indexed::printtree(std::ostream & os, unsigned indent) const { - debugmsg("indexed printtree",LOGLEVEL_PRINT); - os << std::string(indent,' ') << "indexed: " << seq.size() << " indices"; + debugmsg("indexed printtree", LOGLEVEL_PRINT); + GINAC_ASSERT(seq.size() > 0); + + os << std::string(indent, ' ') << class_name() << ", " << seq.size()-1 << " indices"; os << ",hash=" << hashvalue << ",flags=" << flags << std::endl; - printtreeindices(os,indent); + printtreeindices(os, indent); } void indexed::print(std::ostream & os, unsigned upper_precedence) const { - debugmsg("indexed print",LOGLEVEL_PRINT); - os << "UNNAMEDINDEX"; + debugmsg("indexed print", LOGLEVEL_PRINT); + GINAC_ASSERT(seq.size() > 0); + + const ex & base = seq[0]; + bool need_parens = is_ex_exactly_of_type(base, add) || is_ex_exactly_of_type(base, mul) + || is_ex_exactly_of_type(base, ncmul) || is_ex_exactly_of_type(base, power); + if (need_parens) + os << "("; + os << base; + if (need_parens) + os << ")"; printindices(os); } -void indexed::printcsrc(std::ostream & os, unsigned type, - unsigned upper_precedence) const -{ - debugmsg("indexed print csrc",LOGLEVEL_PRINT); - print(os,upper_precedence); -} - bool indexed::info(unsigned inf) const { - if (inf==info_flags::indexed) return true; - if (inf==info_flags::has_indices) return seq.size()!=0; + if (inf == info_flags::indexed) return true; + if (inf == info_flags::has_indices) return seq.size() > 1; return inherited::info(inf); } -// protected - -/** Implementation of ex::diff() for an indexed object. It always returns 0. - * @see ex::diff */ -ex indexed::derivative(const symbol & s) const +bool indexed::all_index_values_are(unsigned inf) const { - return _ex0(); + // No indices? Then no property can be fulfilled + if (seq.size() < 2) + return false; + + // Check all indices + exvector::const_iterator it = seq.begin() + 1, itend = seq.end(); + while (it != itend) { + GINAC_ASSERT(is_ex_of_type(*it, idx)); + if (!ex_to_idx(*it).get_value().info(inf)) + return false; + it++; + } + return true; } int indexed::compare_same_type(const basic & other) const { - GINAC_ASSERT(is_of_type(other,indexed)); + GINAC_ASSERT(is_of_type(other, indexed)); return inherited::compare_same_type(other); } -bool indexed::is_equal_same_type(const basic & other) const +// The main difference between sort_index_vector() and canonicalize_indices() +// is that the latter takes the symmetry of the object into account. Once we +// implement mixed symmetries, canonicalize_indices() will only be able to +// reorder index pairs with known symmetry properties, while sort_index_vector() +// always sorts the whole vector. + +/** Bring a vector of indices into a canonic order (don't care about the + * symmetry of the objects carrying the indices). Dummy indices will lie + * next to each other after the sorting. + * + * @param v Index vector to be sorted */ +static void sort_index_vector(exvector &v) { - GINAC_ASSERT(is_of_type(other,indexed)); - return inherited::is_equal_same_type(other); + // Nothing to sort if less than 2 elements + if (v.size() < 2) + return; + + // Simple bubble sort algorithm should be sufficient for the small + // number of indices expected + exvector::iterator it1 = v.begin(), itend = v.end(), next_to_last_idx = itend - 1; + while (it1 != next_to_last_idx) { + exvector::iterator it2 = it1 + 1; + while (it2 != itend) { + if (it1->compare(*it2) > 0) + it1->swap(*it2); + it2++; + } + it1++; + } } -unsigned indexed::return_type(void) const +/** Bring a vector of indices into a canonic order. This operation only makes + * sense if the object carrying these indices is either symmetric or totally + * antisymmetric with respect to the indices. + * + * @param itbegin Start of index vector + * @param itend End of index vector + * @param antisymm Whether the object is antisymmetric + * @return the sign introduced by the reordering of the indices if the object + * is antisymmetric (or 0 if two equal indices are encountered). For + * symmetric objects, this is always +1. If the index vector was + * already in a canonic order this function returns INT_MAX. */ +static int canonicalize_indices(exvector::iterator itbegin, exvector::iterator itend, bool antisymm) { - return return_types::noncommutative; + bool something_changed = false; + int sig = 1; + + // Simple bubble sort algorithm should be sufficient for the small + // number of indices expected + exvector::iterator it1 = itbegin, next_to_last_idx = itend - 1; + while (it1 != next_to_last_idx) { + exvector::iterator it2 = it1 + 1; + while (it2 != itend) { + int cmpval = it1->compare(*it2); + if (cmpval == 1) { + it1->swap(*it2); + something_changed = true; + if (antisymm) + sig = -sig; + } else if (cmpval == 0 && antisymm) { + something_changed = true; + sig = 0; + } + it2++; + } + it1++; + } + + return something_changed ? sig : INT_MAX; } - -unsigned indexed::return_type_tinfo(void) const + +ex indexed::eval(int level) const { - return tinfo_key; + // First evaluate children, then we will end up here again + if (level > 1) + return indexed(symmetry, evalchildren(level)); + + // Canonicalize indices according to the symmetry properties + if (seq.size() > 2 && (symmetry != unknown && symmetry != mixed)) { + exvector v = seq; + int sig = canonicalize_indices(v.begin() + 1, v.end(), symmetry == antisymmetric); + if (sig != INT_MAX) { + // Something has changed while sorting indices, more evaluations later + if (sig == 0) + return _ex0(); + return ex(sig) * thisexprseq(v); + } + } + + // Let the class of the base object perform additional evaluations + return op(0).bp->eval_indexed(*this); } ex indexed::thisexprseq(const exvector & v) const { - return indexed(v); + return indexed(symmetry, v); } ex indexed::thisexprseq(exvector * vp) const { - return indexed(vp); + return indexed(symmetry, vp); +} + +ex indexed::expand(unsigned options) const +{ + GINAC_ASSERT(seq.size() > 0); + + if ((options & expand_options::expand_indexed) && is_ex_exactly_of_type(seq[0], add)) { + + // expand_indexed expands (a+b).i -> a.i + b.i + const ex & base = seq[0]; + ex sum = _ex0(); + for (unsigned i=0; i 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 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 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; isecond; +} + +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 diff --git a/ginac/indexed.h b/ginac/indexed.h index 94ba0e68..b15ed8b3 100644 --- a/ginac/indexed.h +++ b/ginac/indexed.h @@ -1,6 +1,6 @@ /** @file indexed.h * - * Interface to GiNaC's index carrying objects. */ + * Interface to GiNaC's indexed expressions. */ /* * GiNaC Copyright (C) 1999-2001 Johannes Gutenberg University Mainz, Germany @@ -23,65 +23,229 @@ #ifndef __GINAC_INDEXED_H__ #define __GINAC_INDEXED_H__ -#include +#include + #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 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 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(*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__ diff --git a/ginac/matrix.cpp b/ginac/matrix.cpp index 0c32a767..a63275f5 100644 --- a/ginac/matrix.cpp +++ b/ginac/matrix.cpp @@ -28,6 +28,8 @@ #include "archive.h" #include "numeric.h" #include "lst.h" +#include "idx.h" +#include "indexed.h" #include "utils.h" #include "debugmsg.h" #include "power.h" @@ -157,7 +159,7 @@ void matrix::print(std::ostream & os, unsigned upper_precedence) const void matrix::printraw(std::ostream & os) const { debugmsg("matrix printraw",LOGLEVEL_PRINT); - os << "matrix(" << row << "," << col <<","; + os << class_name() << "(" << row << "," << col <<","; for (unsigned r=0; r(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 ////////// diff --git a/ginac/matrix.h b/ginac/matrix.h index 02203e05..7efbfed7 100644 --- a/ginac/matrix.h +++ b/ginac/matrix.h @@ -50,7 +50,9 @@ public: bool has(const ex & other) const; ex eval(int level=0) const; ex evalf(int level=0) const; - // ex subs(const lst & ls, const lst & lr) const; + ex subs(const lst & ls, const lst & lr) const; + ex eval_indexed(const basic & i) const; + bool contract_with(ex & self, ex & other) const; protected: unsigned return_type(void) const { return return_types::noncommutative; }; // new virtual functions which can be overridden by derived classes diff --git a/ginac/mul.cpp b/ginac/mul.cpp index 887a9169..de5d03c0 100644 --- a/ginac/mul.cpp +++ b/ginac/mul.cpp @@ -418,19 +418,6 @@ ex mul::evalf(int level) const return mul(s,overall_coeff.evalf(level)); } -exvector mul::get_indices(void) const -{ - // return union of indices of factors - exvector iv; - for (epvector::const_iterator cit=seq.begin(); cit!=seq.end(); ++cit) { - exvector subiv=(*cit).rest.get_indices(); - iv.reserve(iv.size()+subiv.size()); - for (exvector::const_iterator cit2=subiv.begin(); cit2!=subiv.end(); ++cit2) - iv.push_back(*cit2); - } - return iv; -} - ex mul::simplify_ncmul(const exvector & v) const { throw(std::logic_error("mul::simplify_ncmul() should never have been called!")); diff --git a/ginac/mul.h b/ginac/mul.h index dba5daf1..c531fd6f 100644 --- a/ginac/mul.h +++ b/ginac/mul.h @@ -61,7 +61,7 @@ public: numeric integer_content(void) const; ex smod(const numeric &xi) const; numeric max_coefficient(void) const; - exvector get_indices(void) const; + exvector get_free_indices(void) const; ex simplify_ncmul(const exvector & v) const; protected: ex derivative(const symbol & s) const; diff --git a/ginac/ncmul.cpp b/ginac/ncmul.cpp index 745cf96b..1433a3a8 100644 --- a/ginac/ncmul.cpp +++ b/ginac/ncmul.cpp @@ -471,20 +471,6 @@ ex ncmul::eval(int level) const status_flags::evaluated); } -exvector ncmul::get_indices(void) const -{ - // return union of indices of factors - exvector iv; - for (exvector::const_iterator cit=seq.begin(); cit!=seq.end(); ++cit) { - exvector subiv=(*cit).get_indices(); - iv.reserve(iv.size()+subiv.size()); - for (exvector::const_iterator cit2=subiv.begin(); cit2!=subiv.end(); ++cit2) { - iv.push_back(*cit2); - } - } - return iv; -} - ex ncmul::subs(const lst & ls, const lst & lr) const { return ncmul(subschildren(ls, lr)); diff --git a/ginac/ncmul.h b/ginac/ncmul.h index cc186a0a..dd96a891 100644 --- a/ginac/ncmul.h +++ b/ginac/ncmul.h @@ -63,7 +63,7 @@ public: ex coeff(const symbol & s, int n=1) const; ex eval(int level=0) const; ex subs(const lst & ls, const lst & lr) const; - exvector get_indices(void) const; + exvector get_free_indices(void) const; ex thisexprseq(const exvector & v) const; ex thisexprseq(exvector * vp) const; protected: diff --git a/ginac/numeric.cpp b/ginac/numeric.cpp index b40e71fc..3bd24a6a 100644 --- a/ginac/numeric.cpp +++ b/ginac/numeric.cpp @@ -449,7 +449,7 @@ void numeric::printraw(std::ostream &os) const // The method printraw doesn't do much, it simply uses CLN's operator<<() // for output, which is ugly but reliable. e.g: 2+2i debugmsg("numeric printraw", LOGLEVEL_PRINT); - os << "numeric(" << cln::the(value) << ")"; + os << class_name() << "(" << cln::the(value) << ")"; } diff --git a/ginac/power.cpp b/ginac/power.cpp index 5ebeae63..b0f9cd2d 100644 --- a/ginac/power.cpp +++ b/ginac/power.cpp @@ -135,7 +135,7 @@ void power::printraw(std::ostream & os) const { debugmsg("power printraw",LOGLEVEL_PRINT); - os << "power("; + os << class_name() << "("; basis.printraw(os); os << ","; exponent.printraw(os); @@ -146,8 +146,8 @@ void power::printtree(std::ostream & os, unsigned indent) const { debugmsg("power printtree",LOGLEVEL_PRINT); - os << std::string(indent,' ') << "power: " - << "hash=" << hashvalue + os << std::string(indent,' ') << class_name() + << ", hash=" << hashvalue << " (0x" << std::hex << hashvalue << std::dec << ")" << ", flags=" << flags << std::endl; basis.printtree(os, indent+delta_indent); diff --git a/ginac/power.h b/ginac/power.h index 8a9d65f8..2c748616 100644 --- a/ginac/power.h +++ b/ginac/power.h @@ -64,6 +64,7 @@ public: ex subs(const lst & ls, const lst & lr) const; ex normal(lst &sym_lst, lst &repl_lst, int level = 0) const; ex to_rational(lst &repl_lst) const; + exvector get_free_indices(void) const; ex simplify_ncmul(const exvector & v) const; protected: ex derivative(const symbol & s) const; diff --git a/ginac/pseries.cpp b/ginac/pseries.cpp index 692efbcd..3b2ac701 100644 --- a/ginac/pseries.cpp +++ b/ginac/pseries.cpp @@ -175,7 +175,7 @@ void pseries::print(std::ostream &os, unsigned upper_precedence) const void pseries::printraw(std::ostream &os) const { debugmsg("pseries printraw", LOGLEVEL_PRINT); - os << "pseries(" << var << ";" << point << ";"; + os << class_name() << "(" << var << ";" << point << ";"; for (epvector::const_iterator i=seq.begin(); i!=seq.end(); ++i) os << "(" << (*i).rest << "," << (*i).coeff << "),"; os << ")"; @@ -185,7 +185,7 @@ void pseries::printraw(std::ostream &os) const void pseries::printtree(std::ostream & os, unsigned indent) const { debugmsg("pseries printtree",LOGLEVEL_PRINT); - os << std::string(indent,' ') << "pseries " + os << std::string(indent,' ') << class_name() << ", hash=" << hashvalue << " (0x" << std::hex << hashvalue << std::dec << ")" << ", flags=" << flags << std::endl; @@ -227,6 +227,7 @@ int pseries::compare_same_type(const basic & other) const ++it; ++o_it; } + // so they are equal. return 0; } diff --git a/ginac/relational.cpp b/ginac/relational.cpp index b22a6387..20ff40b8 100644 --- a/ginac/relational.cpp +++ b/ginac/relational.cpp @@ -143,7 +143,7 @@ void relational::print(std::ostream & os, unsigned upper_precedence) const void relational::printraw(std::ostream & os) const { debugmsg("relational printraw",LOGLEVEL_PRINT); - os << "RELATIONAL("; + os << class_name() << "("; lh.printraw(os); os << ","; rh.printraw(os); diff --git a/ginac/structure.cpp b/ginac/structure.cpp index 35b53e3c..e9c3b021 100644 --- a/ginac/structure.cpp +++ b/ginac/structure.cpp @@ -92,14 +92,14 @@ void structure::printraw(std::ostream & os) const { debugmsg("structure printraw",LOGLEVEL_PRINT); - os << "structure(hash=" << hashvalue << ",flags=" << flags << ")"; + os << class_name() << "(hash=" << hashvalue << ",flags=" << flags << ")"; } void structure::print(std::ostream & os, unsigned upper_precedence) const { debugmsg("structure print",LOGLEVEL_PRINT); - os << "structure()"; + os << class_name() << "()"; } void structure::printtree(std::ostream & os, unsigned indent) const diff --git a/ginac/structure.pl b/ginac/structure.pl index e46bf4b2..9a734c90 100755 --- a/ginac/structure.pl +++ b/ginac/structure.pl @@ -331,13 +331,13 @@ basic * ${STRUCTURE}::duplicate() const void ${STRUCTURE}::printraw(ostream & os) const { debugmsg("${STRUCTURE} printraw",LOGLEVEL_PRINT); - os << "${STRUCTURE}()"; + os << class_name() << "()"; } void ${STRUCTURE}::print(ostream & os, unsigned upper_precedence) const { debugmsg("${STRUCTURE} print",LOGLEVEL_PRINT); - os << "${STRUCTURE}()"; + os << class_name() << "()"; } void ${STRUCTURE}::printtree(ostream & os, unsigned indent) const diff --git a/ginac/symbol.cpp b/ginac/symbol.cpp index 169ee631..4328d19b 100644 --- a/ginac/symbol.cpp +++ b/ginac/symbol.cpp @@ -25,7 +25,6 @@ #include "symbol.h" #include "lst.h" -#include "idx.h" #include "archive.h" #include "debugmsg.h" #include "utils.h" @@ -141,7 +140,7 @@ void symbol::print(std::ostream & os, unsigned upper_precedence) const void symbol::printraw(std::ostream & os) const { debugmsg("symbol printraw",LOGLEVEL_PRINT); - os << "symbol(" << "name=" << name << ",serial=" << serial + os << class_name() << "(" << "name=" << name << ",serial=" << serial << ",hash=" << hashvalue << ",flags=" << flags << ")"; } @@ -225,11 +224,6 @@ ex symbol::eval(int level) const ex symbol::subs(const lst & ls, const lst & lr) const { GINAC_ASSERT(ls.nops()==lr.nops()); -#ifdef DO_GINAC_ASSERT - for (unsigned i=0; i + +#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(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(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; iis_symbolic()) { + for (int i=1; i(*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__ diff --git a/ginac/tinfos.h b/ginac/tinfos.h index 927cf79c..2b189595 100644 --- a/ginac/tinfos.h +++ b/ginac/tinfos.h @@ -25,48 +25,47 @@ 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 -- 2.45.0