#include "mul.h"
#include "ncmul.h"
#include "power.h"
+#include "symmetry.h"
#include "lst.h"
#include "print.h"
#include "archive.h"
// default constructor, destructor, copy constructor assignment operator and helpers
//////////
-indexed::indexed() : symmetry(unknown)
+indexed::indexed() : symtree(sy_none())
{
debugmsg("indexed default constructor", LOGLEVEL_CONSTRUCT);
tinfo_key = TINFO_indexed;
void indexed::copy(const indexed & other)
{
inherited::copy(other);
- symmetry = other.symmetry;
+ symtree = other.symtree;
}
DEFAULT_DESTROY(indexed)
// other constructors
//////////
-indexed::indexed(const ex & b) : inherited(b), symmetry(unknown)
+indexed::indexed(const ex & b) : inherited(b), symtree(sy_none())
{
debugmsg("indexed constructor from ex", LOGLEVEL_CONSTRUCT);
tinfo_key = TINFO_indexed;
- assert_all_indices_of_type_idx();
+ validate();
}
-indexed::indexed(const ex & b, const ex & i1) : inherited(b, i1), symmetry(unknown)
+indexed::indexed(const ex & b, const ex & i1) : inherited(b, i1), symtree(sy_none())
{
debugmsg("indexed constructor from ex,ex", LOGLEVEL_CONSTRUCT);
tinfo_key = TINFO_indexed;
- assert_all_indices_of_type_idx();
+ validate();
}
-indexed::indexed(const ex & b, const ex & i1, const ex & i2) : inherited(b, i1, i2), symmetry(unknown)
+indexed::indexed(const ex & b, const ex & i1, const ex & i2) : inherited(b, i1, i2), symtree(sy_none())
{
debugmsg("indexed constructor from ex,ex,ex", LOGLEVEL_CONSTRUCT);
tinfo_key = TINFO_indexed;
- assert_all_indices_of_type_idx();
+ validate();
}
-indexed::indexed(const ex & b, const ex & i1, const ex & i2, const ex & i3) : inherited(b, i1, i2, i3), symmetry(unknown)
+indexed::indexed(const ex & b, const ex & i1, const ex & i2, const ex & i3) : inherited(b, i1, i2, i3), symtree(sy_none())
{
debugmsg("indexed constructor from ex,ex,ex,ex", LOGLEVEL_CONSTRUCT);
tinfo_key = TINFO_indexed;
- assert_all_indices_of_type_idx();
+ validate();
}
-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)
+indexed::indexed(const ex & b, const ex & i1, const ex & i2, const ex & i3, const ex & i4) : inherited(b, i1, i2, i3, i4), symtree(sy_none())
{
debugmsg("indexed constructor from ex,ex,ex,ex,ex", LOGLEVEL_CONSTRUCT);
tinfo_key = TINFO_indexed;
- assert_all_indices_of_type_idx();
+ validate();
}
-indexed::indexed(const ex & b, symmetry_type symm, const ex & i1, const ex & i2) : inherited(b, i1, i2), symmetry(symm)
+indexed::indexed(const ex & b, const symmetry & symm, const ex & i1, const ex & i2) : inherited(b, i1, i2), symtree(symm)
{
debugmsg("indexed constructor from ex,symmetry,ex,ex", LOGLEVEL_CONSTRUCT);
tinfo_key = TINFO_indexed;
- assert_all_indices_of_type_idx();
+ validate();
}
-indexed::indexed(const ex & b, symmetry_type symm, const ex & i1, const ex & i2, const ex & i3) : inherited(b, i1, i2, i3), symmetry(symm)
+indexed::indexed(const ex & b, const symmetry & symm, const ex & i1, const ex & i2, const ex & i3) : inherited(b, i1, i2, i3), symtree(symm)
{
debugmsg("indexed constructor from ex,symmetry,ex,ex,ex", LOGLEVEL_CONSTRUCT);
tinfo_key = TINFO_indexed;
- assert_all_indices_of_type_idx();
+ validate();
}
-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)
+indexed::indexed(const ex & b, const symmetry & symm, const ex & i1, const ex & i2, const ex & i3, const ex & i4) : inherited(b, i1, i2, i3, i4), symtree(symm)
{
debugmsg("indexed constructor from ex,symmetry,ex,ex,ex,ex", LOGLEVEL_CONSTRUCT);
tinfo_key = TINFO_indexed;
- assert_all_indices_of_type_idx();
+ validate();
}
-indexed::indexed(const ex & b, const exvector & v) : inherited(b), symmetry(unknown)
+indexed::indexed(const ex & b, const exvector & v) : inherited(b), symtree(sy_none())
{
debugmsg("indexed constructor from ex,exvector", LOGLEVEL_CONSTRUCT);
seq.insert(seq.end(), v.begin(), v.end());
tinfo_key = TINFO_indexed;
- assert_all_indices_of_type_idx();
+ validate();
}
-indexed::indexed(const ex & b, symmetry_type symm, const exvector & v) : inherited(b), symmetry(symm)
+indexed::indexed(const ex & b, const symmetry & symm, const exvector & v) : inherited(b), symtree(symm)
{
debugmsg("indexed constructor from ex,symmetry,exvector", LOGLEVEL_CONSTRUCT);
seq.insert(seq.end(), v.begin(), v.end());
tinfo_key = TINFO_indexed;
- assert_all_indices_of_type_idx();
+ validate();
}
-indexed::indexed(symmetry_type symm, const exprseq & es) : inherited(es), symmetry(symm)
+indexed::indexed(const symmetry & symm, const exprseq & es) : inherited(es), symtree(symm)
{
debugmsg("indexed constructor from symmetry,exprseq", LOGLEVEL_CONSTRUCT);
tinfo_key = TINFO_indexed;
- assert_all_indices_of_type_idx();
}
-indexed::indexed(symmetry_type symm, const exvector & v, bool discardable) : inherited(v, discardable), symmetry(symm)
+indexed::indexed(const symmetry & symm, const exvector & v, bool discardable) : inherited(v, discardable), symtree(symm)
{
debugmsg("indexed constructor from symmetry,exvector", LOGLEVEL_CONSTRUCT);
tinfo_key = TINFO_indexed;
- assert_all_indices_of_type_idx();
}
-indexed::indexed(symmetry_type symm, exvector * vp) : inherited(vp), symmetry(symm)
+indexed::indexed(const symmetry & symm, exvector * vp) : inherited(vp), symtree(symm)
{
debugmsg("indexed constructor from symmetry,exvector *", LOGLEVEL_CONSTRUCT);
tinfo_key = TINFO_indexed;
- assert_all_indices_of_type_idx();
}
//////////
indexed::indexed(const archive_node &n, const lst &sym_lst) : inherited(n, sym_lst)
{
debugmsg("indexed constructor from archive_node", LOGLEVEL_CONSTRUCT);
- unsigned int symm;
- if (!(n.find_unsigned("symmetry", symm)))
- throw (std::runtime_error("unknown indexed symmetry type in archive"));
+ if (!n.find_ex("symmetry", symtree, sym_lst)) {
+ // GiNaC versions <= 0.9.0 had an unsigned "symmetry" property
+ unsigned symm = 0;
+ n.find_unsigned("symmetry", symm);
+ switch (symm) {
+ case 1:
+ symtree = sy_symm();
+ break;
+ case 2:
+ symtree = sy_anti();
+ break;
+ default:
+ symtree = sy_none();
+ break;
+ }
+ ex_to_nonconst_symmetry(symtree).validate(seq.size() - 1);
+ }
}
void indexed::archive(archive_node &n) const
{
inherited::archive(n);
- n.add_unsigned("symmetry", symmetry);
+ n.add_ex("symmetry", symtree);
}
DEFAULT_UNARCHIVE(indexed)
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;
- }
+ << ", " << seq.size()-1 << " indices"
+ << ", symmetry=" << symtree << std::endl;
c.s << std::endl;
unsigned delta_indent = static_cast<const print_tree &>(c).delta_indent;
seq[0].print(c, level + delta_indent);
return inherited::info(inf);
}
+struct idx_is_not : public std::binary_function<ex, unsigned, bool> {
+ bool operator() (const ex & e, unsigned inf) const {
+ return !(ex_to_idx(e).get_value().info(inf));
+ }
+};
+
bool indexed::all_index_values_are(unsigned inf) const
{
// No indices? Then no property can be fulfilled
return false;
// Check all indices
- exvector::const_iterator it = seq.begin() + 1, itend = seq.end();
- while (it != itend) {
- GINAC_ASSERT(is_ex_of_type(*it, idx));
- if (!ex_to_idx(*it).get_value().info(inf))
- return false;
- it++;
- }
- return true;
+ return find_if(seq.begin() + 1, seq.end(), bind2nd(idx_is_not(), inf)) == seq.end();
}
int indexed::compare_same_type(const basic & other) const
return inherited::compare_same_type(other);
}
-// The main difference between sort_index_vector() and canonicalize_indices()
-// is that the latter takes the symmetry of the object into account. Once we
-// implement mixed symmetries, canonicalize_indices() will only be able to
-// reorder index pairs with known symmetry properties, while sort_index_vector()
-// always sorts the whole vector.
-
-/** Bring a vector of indices into a canonic order. This operation only makes
- * sense if the object carrying these indices is either symmetric or totally
- * antisymmetric with respect to the indices.
- *
- * @param itbegin Start of index vector
- * @param itend End of index vector
- * @param antisymm Whether the object is antisymmetric
- * @return the sign introduced by the reordering of the indices if the object
- * is antisymmetric (or 0 if two equal indices are encountered). For
- * symmetric objects, this is always +1. If the index vector was
- * already in a canonic order this function returns INT_MAX. */
-static int canonicalize_indices(exvector::iterator itbegin, exvector::iterator itend, bool antisymm)
-{
- bool something_changed = false;
- int sig = 1;
-
- // Simple bubble sort algorithm should be sufficient for the small
- // number of indices expected
- exvector::iterator it1 = itbegin, next_to_last_idx = itend - 1;
- while (it1 != next_to_last_idx) {
- exvector::iterator it2 = it1 + 1;
- while (it2 != itend) {
- int cmpval = it1->compare(*it2);
- if (cmpval == 1) {
- it1->swap(*it2);
- something_changed = true;
- if (antisymm)
- sig = -sig;
- } else if (cmpval == 0 && antisymm) {
- something_changed = true;
- sig = 0;
- }
- it2++;
- }
- it1++;
- }
-
- return something_changed ? sig : INT_MAX;
-}
-
ex indexed::eval(int level) const
{
// First evaluate children, then we will end up here again
if (level > 1)
- return indexed(symmetry, evalchildren(level));
+ return indexed(ex_to_symmetry(symtree), evalchildren(level));
const ex &base = seq[0];
// 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;
+ 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)) {
+ if (seq.size() > 2) {
exvector v = seq;
- int sig = canonicalize_indices(v.begin() + 1, v.end(), symmetry == antisymmetric);
+ GINAC_ASSERT(is_ex_exactly_of_type(symtree, symmetry));
+ int sig = canonicalize(v.begin() + 1, ex_to_symmetry(symtree));
if (sig != INT_MAX) {
// Something has changed while sorting indices, more evaluations later
if (sig == 0)
ex indexed::thisexprseq(const exvector & v) const
{
- return indexed(symmetry, v);
+ return indexed(ex_to_symmetry(symtree), v);
}
ex indexed::thisexprseq(exvector * vp) const
{
- return indexed(symmetry, vp);
+ return indexed(ex_to_symmetry(symtree), vp);
}
ex indexed::expand(unsigned options) 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. */
-void indexed::assert_all_indices_of_type_idx(void) const
+/** Check whether all indices are of class idx and validate the symmetry
+ * tree. This function is used internally to make sure that all constructed
+ * indexed objects really carry indices and not some other classes. */
+void indexed::validate(void) const
{
GINAC_ASSERT(seq.size() > 0);
exvector::const_iterator it = seq.begin() + 1, itend = seq.end();
throw(std::invalid_argument("indices of indexed object must be of type idx"));
it++;
}
+
+ if (!symtree.is_zero()) {
+ if (!is_ex_exactly_of_type(symtree, symmetry))
+ throw(std::invalid_argument("symmetry of indexed object must be of type symmetry"));
+ ex_to_nonconst_symmetry(symtree).validate(seq.size() - 1);
+ }
}
//////////
if (v1.size() != v2.size())
return false;
- // And also the indices themselves
- exvector::const_iterator ait = v1.begin(), aitend = v1.end(),
- bit = v2.begin(), bitend = v2.end();
- while (ait != aitend) {
- if (!ait->is_equal(*bit))
- return false;
- ait++; bit++;
- }
- return true;
+ return equal(v1.begin(), v1.end(), v2.begin(), ex_is_equal());
}
exvector indexed::get_indices(void) const
return dummy_indices;
}
+bool indexed::has_dummy_index_for(const ex & i) const
+{
+ exvector::const_iterator it = seq.begin() + 1, itend = seq.end();
+ while (it != itend) {
+ if (is_dummy_pair(*it, i))
+ return true;
+ it++;
+ }
+ return false;
+}
+
exvector indexed::get_free_indices(void) const
{
exvector free_indices, dummy_indices;
return basis.get_free_indices();
}
-/* Function object for STL sort() */
-struct ex_is_less {
- bool operator() (const ex &lh, const ex &rh) const
- {
- return lh.compare(rh) < 0;
- }
-};
-
/** Rename dummy indices in an expression.
*
* @param e Expression to be worked on
if (local_size == 0)
return e;
- sort(local_dummy_indices.begin(), local_dummy_indices.end(), ex_is_less());
-
if (global_size < local_size) {
// More local indices than we encountered before, add the new ones
int remaining = local_size - global_size;
exvector::const_iterator it = local_dummy_indices.begin(), itend = local_dummy_indices.end();
while (it != itend && remaining > 0) {
- exvector::const_iterator git = global_dummy_indices.begin(), gitend = global_dummy_indices.end();
- while (git != gitend) {
- if (it->is_equal(*git))
- goto found;
- git++;
+ if (find_if(global_dummy_indices.begin(), global_dummy_indices.end(), bind2nd(ex_is_equal(), *it)) == global_dummy_indices.end()) {
+ global_dummy_indices.push_back(*it);
+ global_size++;
+ remaining--;
}
- global_dummy_indices.push_back(*it);
- global_size++;
- remaining--;
-found: it++;
+ it++;
}
- sort(global_dummy_indices.begin(), global_dummy_indices.end(), ex_is_less());
}
// Replace index symbols in expression
for (unsigned i=0; i<local_size; i++) {
ex loc_sym = local_dummy_indices[i].op(0);
ex glob_sym = global_dummy_indices[i].op(0);
- if (!loc_sym.is_equal(glob_sym))
+ if (!loc_sym.is_equal(glob_sym)) {
all_equal = false;
- local_syms.append(loc_sym);
- global_syms.append(glob_sym);
+ local_syms.append(loc_sym);
+ global_syms.append(glob_sym);
+ }
}
if (all_equal)
return e;
if (!is_ex_of_type(*it1, indexed))
continue;
+ bool first_noncommutative = (it1->return_type() != return_types::commutative);
+
// Indexed factor found, get free indices and look for contraction
// candidates
exvector free1, dummy1;
if (!is_ex_of_type(*it2, indexed))
continue;
+ bool second_noncommutative = (it2->return_type() != return_types::commutative);
+
// Find free indices of second factor and merge them with free
// indices of first factor
exvector un;
}
// 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();
+ if (dummy.size() > 1
+ && ex_to_symmetry(ex_to_indexed(*it1).symtree).has_symmetry()
+ && ex_to_symmetry(ex_to_indexed(*it2).symtree).has_symmetry()) {
+
+ // Check all pairs of dummy indices
+ for (unsigned idx1=0; idx1<dummy.size()-1; idx1++) {
+ for (unsigned idx2=idx1+1; idx2<dummy.size(); idx2++) {
+
+ // Try and swap the index pair and check whether the
+ // relative sign changed
+ lst subs_lst(dummy[idx1].op(0), dummy[idx2].op(0)), repl_lst(dummy[idx2].op(0), dummy[idx1].op(0));
+ ex swapped1 = it1->subs(subs_lst, repl_lst);
+ ex swapped2 = it2->subs(subs_lst, repl_lst);
+ if (it1->is_equal(swapped1) && it2->is_equal(-swapped2)
+ || it1->is_equal(-swapped1) && it2->is_equal(swapped2)) {
+ free_indices.clear();
+ return _ex0();
+ }
+ }
+ }
}
// Try to contract the first one with the second one
}
if (contracted) {
contraction_done:
- if (non_commutative
+ if (first_noncommutative || second_noncommutative
|| 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)) {
// 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)));
+ ex r = (non_commutative ? ex(ncmul(v, true)) : ex(mul(v)));
return simplify_indexed(r, free_indices, dummy_indices, sp);
}
}
// Find free indices (concatenate them all and call find_free_and_dummy())
- exvector un, local_dummy_indices;
+ // and all dummy indices that appear
+ exvector un, individual_dummy_indices;
it1 = v.begin(); itend = v.end();
while (it1 != itend) {
- exvector free_indices_of_factor = it1->get_free_indices();
+ exvector free_indices_of_factor;
+ if (is_ex_of_type(*it1, indexed)) {
+ exvector dummy_indices_of_factor;
+ find_free_and_dummy(ex_to_indexed(*it1).seq.begin() + 1, ex_to_indexed(*it1).seq.end(), free_indices_of_factor, dummy_indices_of_factor);
+ individual_dummy_indices.insert(individual_dummy_indices.end(), dummy_indices_of_factor.begin(), dummy_indices_of_factor.end());
+ } else
+ free_indices_of_factor = it1->get_free_indices();
un.insert(un.end(), free_indices_of_factor.begin(), free_indices_of_factor.end());
it1++;
}
+ exvector local_dummy_indices;
find_free_and_dummy(un, free_indices, local_dummy_indices);
+ local_dummy_indices.insert(local_dummy_indices.end(), individual_dummy_indices.begin(), individual_dummy_indices.end());
ex r;
if (something_changed)
- r = non_commutative ? ex(ncmul(v)) : ex(mul(v));
+ r = non_commutative ? ex(ncmul(v, true)) : ex(mul(v));
else
r = e;
ex e_expanded = e.expand();
// Simplification of single indexed object: just find the free indices
- // (and perform dummy index renaming if
+ // and perform dummy index renaming
if (is_ex_of_type(e_expanded, indexed)) {
const indexed &i = ex_to_indexed(e_expanded);
exvector local_dummy_indices;
return e_expanded;
}
-ex simplify_indexed(const ex & e)
+/** Simplify/canonicalize expression containing indexed objects. This
+ * performs contraction of dummy indices where possible and checks whether
+ * the free indices in sums are consistent.
+ *
+ * @return simplified expression */
+ex ex::simplify_indexed(void) const
{
exvector free_indices, dummy_indices;
scalar_products sp;
- return simplify_indexed(e, free_indices, dummy_indices, sp);
+ return GiNaC::simplify_indexed(*this, free_indices, dummy_indices, sp);
}
-ex simplify_indexed(const ex & e, const scalar_products & sp)
+/** Simplify/canonicalize expression containing indexed objects. This
+ * performs contraction of dummy indices where possible, checks whether
+ * the free indices in sums are consistent, and automatically replaces
+ * scalar products by known values if desired.
+ *
+ * @param sp Scalar products to be replaced automatically
+ * @return simplified expression */
+ex ex::simplify_indexed(const scalar_products & sp) const
{
exvector free_indices, dummy_indices;
- return simplify_indexed(e, free_indices, dummy_indices, sp);
+ return GiNaC::simplify_indexed(*this, free_indices, dummy_indices, sp);
+}
+
+/** Symmetrize expression over its free indices. */
+ex ex::symmetrize(void) const
+{
+ return GiNaC::symmetrize(*this, get_free_indices());
+}
+
+/** Antisymmetrize expression over its free indices. */
+ex ex::antisymmetrize(void) const
+{
+ return GiNaC::antisymmetrize(*this, get_free_indices());
+}
+
+/** Symmetrize expression by cyclic permutation over its free indices. */
+ex ex::symmetrize_cyclic(void) const
+{
+ return GiNaC::symmetrize_cyclic(*this, get_free_indices());
}
//////////
git://www.ginac.de / ginac.git / blobdiff