symmetry = other.symmetry;
}
-void indexed::destroy(bool call_parent)
-{
- if (call_parent)
- inherited::destroy(call_parent);
-}
+DEFAULT_DESTROY(indexed)
//////////
// other constructors
{
debugmsg("indexed constructor from ex", LOGLEVEL_CONSTRUCT);
tinfo_key = TINFO_indexed;
- GINAC_ASSERT(all_indices_of_type_idx());
+ assert_all_indices_of_type_idx();
}
indexed::indexed(const ex & b, const ex & i1) : inherited(b, i1), symmetry(unknown)
{
debugmsg("indexed constructor from ex,ex", LOGLEVEL_CONSTRUCT);
tinfo_key = TINFO_indexed;
- GINAC_ASSERT(all_indices_of_type_idx());
+ assert_all_indices_of_type_idx();
}
indexed::indexed(const ex & b, const ex & i1, const ex & i2) : inherited(b, i1, i2), symmetry(unknown)
{
debugmsg("indexed constructor from ex,ex,ex", LOGLEVEL_CONSTRUCT);
tinfo_key = TINFO_indexed;
- GINAC_ASSERT(all_indices_of_type_idx());
+ assert_all_indices_of_type_idx();
}
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,ex", LOGLEVEL_CONSTRUCT);
tinfo_key = TINFO_indexed;
- GINAC_ASSERT(all_indices_of_type_idx());
+ assert_all_indices_of_type_idx();
}
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());
+ 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,symmetry,ex,ex", LOGLEVEL_CONSTRUCT);
tinfo_key = TINFO_indexed;
- GINAC_ASSERT(all_indices_of_type_idx());
+ assert_all_indices_of_type_idx();
}
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 ex,symmetry,ex,ex,ex", LOGLEVEL_CONSTRUCT);
tinfo_key = TINFO_indexed;
- GINAC_ASSERT(all_indices_of_type_idx());
+ 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());
+ 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());
+ 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());
+ 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());
+ assert_all_indices_of_type_idx();
}
indexed::indexed(symmetry_type symm, const exvector & v, bool discardable) : inherited(v, discardable), symmetry(symm)
{
debugmsg("indexed constructor from symmetry,exvector", LOGLEVEL_CONSTRUCT);
tinfo_key = TINFO_indexed;
- GINAC_ASSERT(all_indices_of_type_idx());
+ 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());
+ assert_all_indices_of_type_idx();
}
//////////
// 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
//////////
if (level > 1)
return indexed(symmetry, evalchildren(level));
+ const ex &base = seq[0];
+
+ // If the base object is 0, the whole object is 0
+ if (base.is_zero())
+ return _ex0();
+
+ // If the base object is a product, pull out the numeric factor
+ if (is_ex_exactly_of_type(base, mul) && is_ex_exactly_of_type(base.op(base.nops() - 1), numeric)) {
+ exvector v = seq;
+ ex f = ex_to_numeric(base.op(base.nops() - 1));
+ v[0] = seq[0] / f;
+ return f * thisexprseq(v);
+ }
+
// Canonicalize indices according to the symmetry properties
if (seq.size() > 2 && (symmetry != unknown && symmetry != mixed)) {
exvector v = seq;
}
// Let the class of the base object perform additional evaluations
- return op(0).bp->eval_indexed(*this);
+ return base.bp->eval_indexed(*this);
}
ex indexed::thisexprseq(const exvector & v) const
/** 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
+void indexed::assert_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;
+ throw(std::invalid_argument("indices of indexed object must be of type idx"));
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 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;
+ find_free_and_dummy(seq.begin() + 1, seq.end(), 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;
}
}
+ // Contraction of symmetric with antisymmetric object is zero
+ if ((ex_to_indexed(*it1).symmetry == indexed::symmetric &&
+ ex_to_indexed(*it2).symmetry == indexed::antisymmetric
+ || ex_to_indexed(*it1).symmetry == indexed::antisymmetric &&
+ ex_to_indexed(*it2).symmetry == indexed::symmetric)
+ && dummy.size() > 1) {
+ free_indices.clear();
+ return _ex0();
+ }
+
// Try to contract the first one with the second one
- bool contracted = it1->op(0).bp->contract_with(*it1, *it2);
+ bool contracted = it1->op(0).bp->contract_with(it1, it2, v);
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);
+ contracted = it2->op(0).bp->contract_with(it2, it1, v);
}
if (contracted) {
something_changed = true;
}
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)
- return ncmul(v);
+ r = ncmul(v);
else
- return mul(v);
+ r = mul(v);
} else
- return e;
+ r = e;
+
+ // Product of indexed object with a scalar?
+ if (is_ex_exactly_of_type(r, mul) && r.nops() == 2
+ && is_ex_exactly_of_type(r.op(1), numeric) && is_ex_of_type(r.op(0), indexed))
+ return r.op(0).op(0).bp->scalar_mul_indexed(r.op(0), ex_to_numeric(r.op(1)));
+ else
+ return r;
}
/** Simplify indexed expression, return list of free indices. */
// Simplification of sum = sum of simplifications, check consistency of
// free indices in each term
if (is_ex_exactly_of_type(e_expanded, add)) {
+ bool first = true;
ex sum = _ex0();
+ free_indices.clear();
for (unsigned i=0; i<e_expanded.nops(); i++) {
exvector free_indices_of_term;
- sum += simplify_indexed(e_expanded.op(i), free_indices_of_term, sp);
- if (i == 0)
- free_indices = free_indices_of_term;
- else if (!indices_consistent(free_indices, free_indices_of_term))
- throw (std::runtime_error("simplify_indexed: inconsistent indices in sum"));
+ ex term = simplify_indexed(e_expanded.op(i), free_indices_of_term, sp);
+ if (!term.is_zero()) {
+ if (first) {
+ free_indices = free_indices_of_term;
+ sum = term;
+ first = false;
+ } else {
+ if (!indices_consistent(free_indices, free_indices_of_term))
+ throw (std::runtime_error("simplify_indexed: inconsistent indices in sum"));
+ if (is_ex_of_type(sum, indexed) && is_ex_of_type(term, indexed))
+ sum = sum.op(0).bp->add_indexed(sum, term);
+ else
+ sum += term;
+ }
+ }
}
return sum;
git://www.ginac.de / ginac.git / blobdiff