From: Christian Bauer Date: Wed, 9 Oct 2002 20:12:16 +0000 (+0000) Subject: scalar_products allows the specification of an index dimension, e.g. to X-Git-Tag: release_1-1-0~66 X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?a=commitdiff_plain;h=81e5c2c3f18302ef44cff390b790420ce616b83f;p=ginac.git scalar_products allows the specification of an index dimension, e.g. to give different values to 4- and D-dimensional products --- diff --git a/ginac/idx.cpp b/ginac/idx.cpp index 59c3ebc3..d663be71 100644 --- a/ginac/idx.cpp +++ b/ginac/idx.cpp @@ -451,12 +451,7 @@ ex idx::replace_dim(const ex & new_dim) const ex idx::minimal_dim(const idx & other) const { - if (dim.is_equal(other.dim) || dim < other.dim || (is_a(dim) && is_a(other.dim))) - return dim; - else if (dim > other.dim || (is_a(dim) && is_a(other.dim))) - return other.dim; - else - throw (std::runtime_error("idx::minimal_dim: index dimensions cannot be ordered")); + return GiNaC::minimal_dim(dim, other.dim); } ex varidx::toggle_variance(void) const @@ -547,4 +542,14 @@ void find_free_and_dummy(exvector::const_iterator it, exvector::const_iterator i out_free.push_back(*last); } +ex minimal_dim(const ex & dim1, const ex & dim2) +{ + if (dim1.is_equal(dim2) || dim1 < dim2 || (is_a(dim1) && is_a(dim2))) + return dim1; + else if (dim1 > dim2 || (is_a(dim1) && is_a(dim2))) + return dim2; + else + throw (std::runtime_error("minimal_dim(): index dimensions cannot be ordered")); +} + } // namespace GiNaC diff --git a/ginac/idx.h b/ginac/idx.h index 0b31696c..b05d555b 100644 --- a/ginac/idx.h +++ b/ginac/idx.h @@ -263,6 +263,11 @@ inline unsigned count_free_indices(const exvector & v) return free_indices.size(); } +/** Return the minimum of two index dimensions. If this is undecidable, + * throw an exception. Numeric dimensions are always considered "smaller" + * than symbolic dimensions. */ +ex minimal_dim(const ex & dim1, const ex & dim2); + } // namespace GiNaC #endif // ndef __GINAC_IDX_H__ diff --git a/ginac/indexed.cpp b/ginac/indexed.cpp index d1e9efaa..bed59fcd 100644 --- a/ginac/indexed.cpp +++ b/ginac/indexed.cpp @@ -701,8 +701,25 @@ try_again: // At least one dummy index, is it a defined scalar product? bool contracted = false; if (free.empty()) { - if (sp.is_defined(*it1, *it2)) { - *it1 = sp.evaluate(*it1, *it2); + + // Find minimal dimension of all indices of both factors + exvector::const_iterator dit = ex_to(*it1).seq.begin() + 1, ditend = ex_to(*it1).seq.end(); + ex dim = ex_to(*dit).get_dim(); + ++dit; + for (; dit != ditend; ++dit) { + dim = minimal_dim(dim, ex_to(*dit).get_dim()); + } + dit = ex_to(*it2).seq.begin() + 1; + ditend = ex_to(*it2).seq.end(); + for (; dit != ditend; ++dit) { + dim = minimal_dim(dim, ex_to(*dit).get_dim()); + } + + // User-defined scalar product? + if (sp.is_defined(*it1, *it2, dim)) { + + // Yes, substitute it + *it1 = sp.evaluate(*it1, *it2, dim); *it2 = _ex1; goto contraction_done; } @@ -951,12 +968,76 @@ ex ex::symmetrize_cyclic(void) const // helper classes ////////// +spmapkey::spmapkey(const ex & v1_, const ex & v2_, const ex & dim_) : dim(dim_) +{ + // If indexed, extract base objects + ex s1 = is_a(v1_) ? v1_.op(0) : v1_; + ex s2 = is_a(v2_) ? v2_.op(0) : v2_; + + // Enforce canonical order in pair + if (s1.compare(s2) > 0) { + v1 = s2; + v2 = s1; + } else { + v1 = s1; + v2 = s2; + } +} + +bool spmapkey::operator==(const spmapkey &other) const +{ + if (!v1.is_equal(other.v1)) + return false; + if (!v2.is_equal(other.v2)) + return false; + if (is_a(dim) || is_a(other.dim)) + return true; + else + return dim.is_equal(other.dim); +} + +bool spmapkey::operator<(const spmapkey &other) const +{ + int cmp = v1.compare(other.v1); + if (cmp) + return cmp < 0; + cmp = v2.compare(other.v2); + if (cmp) + return cmp < 0; + + // Objects are equal, now check dimensions + if (is_a(dim) || is_a(other.dim)) + return false; + else + return dim.compare(other.dim) < 0; +} + +void spmapkey::debugprint(void) const +{ + std::cerr << "(" << v1 << "," << v2 << "," << dim << ")"; +} + +scalar_products::scalar_products(const scalar_products & other) : spm(other.spm) {} + +const scalar_products & scalar_products::operator=(const scalar_products & other) +{ + if (this != &other) { + spm = other.spm; + } + return *this; +} + void scalar_products::add(const ex & v1, const ex & v2, const ex & sp) { - spm[make_key(v1, v2)] = sp; + spm[spmapkey(v1, v2)] = sp; } -void scalar_products::add_vectors(const lst & l) +void scalar_products::add(const ex & v1, const ex & v2, const ex & dim, const ex & sp) +{ + spm[spmapkey(v1, v2, dim)] = sp; +} + +void scalar_products::add_vectors(const lst & l, const ex & dim) { // Add all possible pairs of products unsigned num = l.nops(); @@ -964,7 +1045,7 @@ void scalar_products::add_vectors(const lst & l) ex a = l.op(i); for (unsigned j=0; jsecond; + return spm.find(spmapkey(v1, v2, dim))->second; } void scalar_products::debugprint(void) const @@ -992,24 +1073,11 @@ void scalar_products::debugprint(void) const spmap::const_iterator i = spm.begin(), end = spm.end(); while (i != end) { const spmapkey & k = i->first; - std::cerr << "item key=(" << k.first << "," << k.second; - std::cerr << "), value=" << i->second << std::endl; + std::cerr << "item key="; + k.debugprint(); + std::cerr << ", value=" << i->second << std::endl; ++i; } } -/** Make key from object pair. */ -spmapkey scalar_products::make_key(const ex & v1, const ex & v2) -{ - // If indexed, extract base objects - ex s1 = is_a(v1) ? v1.op(0) : v1; - ex s2 = is_a(v2) ? 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 ad18ddb5..4d232856 100644 --- a/ginac/indexed.h +++ b/ginac/indexed.h @@ -26,6 +26,7 @@ #include #include "exprseq.h" +#include "wildcard.h" namespace GiNaC { @@ -192,17 +193,22 @@ protected: }; -typedef std::pair spmapkey; +class spmapkey { +public: + spmapkey() : dim(wild()) {} + spmapkey(const ex & v1, const ex & v2, const ex & dim = wild()); + ~spmapkey() {} + + bool operator==(const spmapkey &other) const; + bool operator<(const spmapkey &other) const; -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)); - } + void debugprint(void) const; + +protected: + ex v1, v2, dim; }; -typedef std::map spmap; +typedef std::map spmap; /** Helper class for storing information about known scalar products which * are to be automatically replaced by simplify_indexed(). @@ -210,24 +216,30 @@ typedef std::map spmap; * @see simplify_indexed */ class scalar_products { public: + scalar_products() {} + ~scalar_products() {} + scalar_products(const scalar_products & other); + const scalar_products & operator=(const scalar_products & other); + /** Register scalar product pair and its value. */ void add(const ex & v1, const ex & v2, const ex & sp); + /** Register scalar product pair and its value for a specific space dimension. */ + void add(const ex & v1, const ex & v2, const ex & dim, const ex & sp); + /** Register list of vectors. This adds all possible pairs of products * a.i * b.i with the value a*b (note that this is not a scalar vector * product but an ordinary product of scalars). */ - void add_vectors(const lst & l); + void add_vectors(const lst & l, const ex & dim = wild()); /** 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; + bool is_defined(const ex & v1, const ex & v2, const ex & dim) const; + ex evaluate(const ex & v1, const ex & v2, const ex & dim) const; void debugprint(void) const; -private: - static spmapkey make_key(const ex & v1, const ex & v2); - +protected: spmap spm; /*< Map from defined scalar product pairs to their values */ };