#include "mul.h"
#include "ncmul.h"
#include "power.h"
+#include "lst.h"
+#include "print.h"
#include "archive.h"
#include "utils.h"
#include "debugmsg.h"
symmetry = other.symmetry;
}
-void indexed::destroy(bool call_parent)
-{
- if (call_parent)
- inherited::destroy(call_parent);
-}
+DEFAULT_DESTROY(indexed)
//////////
// other constructors
// archiving
//////////
-/** Construct object from archive_node. */
indexed::indexed(const archive_node &n, const lst &sym_lst) : inherited(n, sym_lst)
{
debugmsg("indexed constructor from archive_node", LOGLEVEL_CONSTRUCT);
throw (std::runtime_error("unknown indexed symmetry type in archive"));
}
-/** Unarchive the object. */
-ex indexed::unarchive(const archive_node &n, const lst &sym_lst)
-{
- return (new indexed(n, sym_lst))->setflag(status_flags::dynallocated);
-}
-
-/** Archive the object. */
void indexed::archive(archive_node &n) const
{
inherited::archive(n);
n.add_unsigned("symmetry", symmetry);
}
+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);
+ 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)
+ c.s << "(";
+ base.print(c);
+ if (need_parens)
+ c.s << ")";
+ printindices(c, level);
+ }
}
bool indexed::info(unsigned inf) 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);
// 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->print(c, level);
it++;
}
}
// 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 two sorted index vectors are consistent (i.e. equal). */
static bool indices_consistent(const exvector & v1, const exvector & v2)
{
return true;
}
+exvector indexed::get_indices(void) const
+{
+ GINAC_ASSERT(seq.size() >= 1);
+ return exvector(seq.begin() + 1, seq.end());
+}
+
exvector indexed::get_dummy_indices(void) const
{
exvector free_indices, dummy_indices;
return dummy_indices;
}
+exvector indexed::get_dummy_indices(const indexed & other) const
+{
+ exvector indices = get_free_indices();
+ exvector other_indices = other.get_free_indices();
+ indices.insert(indices.end(), other_indices.begin(), other_indices.end());
+ exvector dummy_indices;
+ find_dummy_indices(indices, dummy_indices);
+ return dummy_indices;
+}
+
exvector indexed::get_free_indices(void) const
{
exvector free_indices, dummy_indices;
// And remove the dummy indices
exvector free_indices, dummy_indices;
- find_free_and_dummy(un.begin(), un.end(), free_indices, dummy_indices);
+ find_free_and_dummy(un, free_indices, dummy_indices);
return free_indices;
}
// And remove the dummy indices
exvector free_indices, dummy_indices;
- find_free_and_dummy(un.begin(), un.end(), free_indices, dummy_indices);
+ find_free_and_dummy(un, free_indices, dummy_indices);
return free_indices;
}
} else if (is_ex_exactly_of_type(f, ncmul)) {
// Noncommutative factor found, split it as well
non_commutative = true; // everything becomes noncommutative, ncmul will sort out the commutative factors later
- for (int j=0; j<f.nops(); i++)
+ for (int j=0; j<f.nops(); j++)
v.push_back(f.op(j));
} else
v.push_back(f);
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.begin(), un.end(), 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;
}
}
}
// 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
contracted = it2->op(0).bp->contract_with(it2, it1, v);
}
if (contracted) {
- something_changed = true;
+contraction_done:
+ if (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)) {
+
+ // One of the factors became a sum or product:
+ // re-expand expression and run again
+ 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;
}
}
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.begin(), un.end(), free_indices, dummy_indices);
+ 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?