]> www.ginac.de Git - ginac.git/blobdiff - ginac/indexed.cpp
- color and clifford classes are quite functional now
[ginac.git] / ginac / indexed.cpp
index 0ab989e801baccc8a94d60a7df7306da1822a78f..3240c8dd7fe86293918c0aabf807cecd91e101b8 100644 (file)
@@ -28,6 +28,8 @@
 #include "mul.h"
 #include "ncmul.h"
 #include "power.h"
+#include "lst.h"
+#include "print.h"
 #include "archive.h"
 #include "utils.h"
 #include "debugmsg.h"
@@ -175,42 +177,44 @@ DEFAULT_UNARCHIVE(indexed)
 // functions overriding virtual functions from bases classes
 //////////
 
-void indexed::printraw(std::ostream & os) const
+void indexed::print(const print_context & c, unsigned level) const
 {
-       debugmsg("indexed printraw", LOGLEVEL_PRINT);
+       debugmsg("indexed print", 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);
-       GINAC_ASSERT(seq.size() > 0);
+       if (is_of_type(c, print_tree)) {
 
-       os << std::string(indent, ' ') << class_name() << ", " << seq.size()-1 << " indices";
-       os << ",hash=" << hashvalue << ",flags=" << flags << std::endl;
-       printtreeindices(os, indent);
-}
+               c.s << std::string(level, ' ') << class_name()
+                   << std::hex << ", hash=0x" << hashvalue << ", flags=0x" << flags << std::dec
+                   << ", " << seq.size()-1 << " indices";
+               switch (symmetry) {
+                       case symmetric: c.s << ", symmetric"; break;
+                       case antisymmetric: c.s << ", antisymmetric"; break;
+                       default: break;
+               }
+               c.s << std::endl;
+               unsigned delta_indent = static_cast<const print_tree &>(c).delta_indent;
+               seq[0].print(c, level + delta_indent);
+               printindices(c, level + delta_indent);
 
-void indexed::print(std::ostream & os, unsigned upper_precedence) const
-{
-       debugmsg("indexed print", LOGLEVEL_PRINT);
-       GINAC_ASSERT(seq.size() > 0);
+       } else {
 
-       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);
+               bool is_tex = is_of_type(c, print_latex);
+               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)
+                               || is_ex_of_type(base, indexed);
+               if (is_tex)
+                       c.s << "{";
+               if (need_parens)
+                       c.s << "(";
+               base.print(c);
+               if (need_parens)
+                       c.s << ")";
+               if (is_tex)
+                       c.s << "}";
+               printindices(c, level);
+       }
 }
 
 bool indexed::info(unsigned inf) const
@@ -350,6 +354,24 @@ ex indexed::eval(int level) const
        return base.bp->eval_indexed(*this);
 }
 
+int indexed::degree(const ex & s) const
+{
+       return is_equal(*s.bp) ? 1 : 0;
+}
+
+int indexed::ldegree(const ex & s) const
+{
+       return is_equal(*s.bp) ? 1 : 0;
+}
+
+ex indexed::coeff(const ex & s, int n) const
+{
+       if (is_equal(*s.bp))
+               return n==1 ? _ex1() : _ex0();
+       else
+               return n==0 ? ex(*this) : _ex0();
+}
+
 ex indexed::thisexprseq(const exvector & v) const
 {
        return indexed(symmetry, v);
@@ -390,39 +412,43 @@ ex indexed::expand(unsigned options) const
 // non-virtual functions in this class
 //////////
 
-void indexed::printrawindices(std::ostream & os) const
+void indexed::printindices(const print_context & c, unsigned level) const
 {
        if (seq.size() > 1) {
-               exvector::const_iterator it=seq.begin() + 1, itend = seq.end();
-               while (it != itend) {
-                       it->printraw(os);
-                       it++;
-                       if (it != itend)
-                               os << ",";
-               }
-       }
-}
 
-void indexed::printtreeindices(std::ostream & os, unsigned indent) const
-{
-       if (seq.size() > 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() > 1) {
-               exvector::const_iterator it=seq.begin() + 1, itend = seq.end();
-               while (it != itend) {
-                       it->print(os);
-                       it++;
+               if (is_of_type(c, print_latex)) {
+
+                       // TeX output: group by variance
+                       bool first = true;
+                       bool covariant = true;
+
+                       while (it != itend) {
+                               bool cur_covariant = (is_ex_of_type(*it, varidx) ? ex_to_varidx(*it).is_covariant() : true);
+                               if (first || cur_covariant != covariant) {
+                                       if (!first)
+                                               c.s << "}";
+                                       covariant = cur_covariant;
+                                       if (covariant)
+                                               c.s << "_{";
+                                       else
+                                               c.s << "^{";
+                               }
+                               it->print(c, level);
+                               c.s << " ";
+                               first = false;
+                               it++;
+                       }
+                       c.s << "}";
+
+               } else {
+
+                       // Ordinary output
+                       while (it != itend) {
+                               it->print(c, level);
+                               it++;
+                       }
                }
        }
 }
@@ -587,28 +613,36 @@ try_again:
                if (!is_ex_of_type(*it1, indexed))
                        continue;
 
-               // Indexed factor found, look for contraction candidates
+               // Indexed factor found, get free indices and look for contraction
+               // candidates
+               exvector free1, dummy1;
+               find_free_and_dummy(ex_to_indexed(*it1).seq.begin() + 1, ex_to_indexed(*it1).seq.end(), free1, dummy1);
+
                exvector::iterator it2;
                for (it2 = it1 + 1; it2 != itend; it2++) {
 
                        if (!is_ex_of_type(*it2, indexed))
                                continue;
 
+                       // Find free indices of second factor and merge them with free
+                       // indices of first factor
+                       exvector un;
+                       find_free_and_dummy(ex_to_indexed(*it2).seq.begin() + 1, ex_to_indexed(*it2).seq.end(), un, dummy1);
+                       un.insert(un.end(), free1.begin(), free1.end());
+
                        // 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, free, dummy);
                        if (dummy.size() == 0)
                                continue;
 
                        // At least one dummy index, is it a defined scalar product?
+                       bool contracted = false;
                        if (free.size() == 0) {
                                if (sp.is_defined(*it1, *it2)) {
                                        *it1 = sp.evaluate(*it1, *it2);
                                        *it2 = _ex1();
-                                       something_changed = true;
-                                       goto try_again;
+                                       goto contraction_done;
                                }
                        }
 
@@ -623,7 +657,7 @@ try_again:
                        }
 
                        // Try to contract the first one with the second one
-                       bool contracted = it1->op(0).bp->contract_with(it1, it2, v);
+                       contracted = it1->op(0).bp->contract_with(it1, it2, v);
                        if (!contracted) {
 
                                // That didn't work; maybe the second object knows how to
@@ -631,11 +665,25 @@ try_again:
                                contracted = it2->op(0).bp->contract_with(it2, it1, v);
                        }
                        if (contracted) {
-                               something_changed = true;
+contraction_done:
+                               if (non_commutative
+                                || is_ex_exactly_of_type(*it1, add) || is_ex_exactly_of_type(*it2, add)
+                                || is_ex_exactly_of_type(*it1, mul) || is_ex_exactly_of_type(*it2, mul)
+                                || is_ex_exactly_of_type(*it1, ncmul) || is_ex_exactly_of_type(*it2, ncmul)) {
+
+                                       // One of the factors became a sum or product:
+                                       // re-expand expression and run again
+                                       // Non-commutative products are always re-expanded to give
+                                       // simplify_ncmul() the chance to re-order and canonicalize
+                                       // the product
+                                       ex r = (non_commutative ? ex(ncmul(v)) : ex(mul(v)));
+                                       return simplify_indexed(r, free_indices, sp);
+                               }
 
                                // Both objects may have new indices now or they might
                                // even not be indexed objects any more, so we have to
                                // start over
+                               something_changed = true;
                                goto try_again;
                        }
                }
@@ -645,21 +693,16 @@ try_again:
        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());
-               }
+               exvector free_indices_of_factor = it1->get_free_indices();
+               un.insert(un.end(), free_indices_of_factor.begin(), free_indices_of_factor.end());
                it1++;
        }
        find_free_and_dummy(un, free_indices, dummy_indices);
 
        ex r;
-       if (something_changed) {
-               if (non_commutative)
-                       r = ncmul(v);
-               else
-                       r = mul(v);
-       } else
+       if (something_changed)
+               r = non_commutative ? ex(ncmul(v)) : ex(mul(v));
+       else
                r = e;
 
        // Product of indexed object with a scalar?
@@ -746,6 +789,19 @@ void scalar_products::add(const ex & v1, const ex & v2, const ex & sp)
        spm[make_key(v1, v2)] = sp;
 }
 
+void scalar_products::add_vectors(const lst & l)
+{
+       // Add all possible pairs of products
+       unsigned num = l.nops();
+       for (unsigned i=0; i<num; i++) {
+               ex a = l.op(i);
+               for (unsigned j=0; j<num; j++) {
+                       ex b = l.op(j);
+                       add(a, b, a*b);
+               }
+       }
+}
+
 void scalar_products::clear(void)
 {
        spm.clear();