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 59c3ebc3eccd67938d8b6d8714201c06fed3a40a..d663be710e7ec0b59dfe89344986e7763adc55c0 100644 (file)
--- 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<numeric>(dim) && is_a<symbol>(other.dim)))
-               return dim;
-       else if (dim > other.dim || (is_a<symbol>(dim) && is_a<numeric>(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<numeric>(dim1) && is_a<symbol>(dim2)))
+               return dim1;
+       else if (dim1 > dim2 || (is_a<symbol>(dim1) && is_a<numeric>(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 0b31696c772840c347aa0fad1aa9e71099de93c7..b05d555b2a26b05f01b8ea90c3edfed6354ed7d4 100644 (file)
--- 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 d1e9efaaeb02fc371f11a810b53046ec07378273..bed59fcda2e0d7cfe7dd5bbd384e1f549130cde2 100644 (file)
--- 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<indexed>(*it1).seq.begin() + 1, ditend = ex_to<indexed>(*it1).seq.end();
+                               ex dim = ex_to<idx>(*dit).get_dim();
+                               ++dit;
+                               for (; dit != ditend; ++dit) {
+                                       dim = minimal_dim(dim, ex_to<idx>(*dit).get_dim());
+                               }
+                               dit = ex_to<indexed>(*it2).seq.begin() + 1;
+                               ditend = ex_to<indexed>(*it2).seq.end();
+                               for (; dit != ditend; ++dit) {
+                                       dim = minimal_dim(dim, ex_to<idx>(*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<indexed>(v1_) ? v1_.op(0) : v1_;
+       ex s2 = is_a<indexed>(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<wildcard>(dim) || is_a<wildcard>(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<wildcard>(dim) || is_a<wildcard>(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; j<num; j++) {
                        ex b = l.op(j);
-                       add(a, b, a*b);
+                       add(a, b, dim, a*b);
                }
        }
 }
@@ -975,15 +1056,15 @@ void scalar_products::clear(void)
 }
 
 /** Check whether scalar product pair is defined. */
-bool scalar_products::is_defined(const ex & v1, const ex & v2) const
+bool scalar_products::is_defined(const ex & v1, const ex & v2, const ex & dim) const
 {
-       return spm.find(make_key(v1, v2)) != spm.end();
+       return spm.find(spmapkey(v1, v2, dim)) != spm.end();
 }
 
 /** Return value of defined scalar product pair. */
-ex scalar_products::evaluate(const ex & v1, const ex & v2) const
+ex scalar_products::evaluate(const ex & v1, const ex & v2, const ex & dim) const
 {
-       return spm.find(make_key(v1, v2))->second;
+       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<indexed>(v1) ? v1.op(0) : v1;
-       ex s2 = is_a<indexed>(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 ad18ddb55f91c568e1b09f727a4f8b814ff59329..4d232856e7b078d160eea325df2e49dfb6b05add 100644 (file)
--- a/ginac/indexed.h
+++ b/ginac/indexed.h
@@ -26,6 +26,7 @@
 #include <map>
 
 #include "exprseq.h"
+#include "wildcard.h"
 
 namespace GiNaC {
 
@@ -192,17 +193,22 @@ protected:
 };
 
 
-typedef std::pair<ex, ex> 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<spmapkey, ex, spmapkey_is_less> spmap;
+typedef std::map<spmapkey, ex> 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<spmapkey, ex, spmapkey_is_less> 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 */
 };