From: Christian Bauer Date: Mon, 11 Jun 2001 23:49:51 +0000 (+0000) Subject: introduced new class for constructing symmetry tree definitions X-Git-Tag: release_0-9-1~43 X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=commitdiff_plain;h=5f1c3b3861f1fc7978b5f734e6e058ba95de355c introduced new class for constructing symmetry tree definitions --- diff --git a/ginac/symmetry.cpp b/ginac/symmetry.cpp new file mode 100644 index 00000000..db055808 --- /dev/null +++ b/ginac/symmetry.cpp @@ -0,0 +1,412 @@ +/** @file symmetry.cpp + * + * Implementation of GiNaC's symmetry definitions. */ + +/* + * GiNaC Copyright (C) 1999-2001 Johannes Gutenberg University Mainz, Germany + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + */ + +#include +#include +#include + +#define DO_GINAC_ASSERT +#include "assertion.h" + +#include "symmetry.h" +#include "lst.h" +#include "numeric.h" // for factorial() +#include "print.h" +#include "archive.h" +#include "utils.h" +#include "debugmsg.h" + +namespace GiNaC { + +GINAC_IMPLEMENT_REGISTERED_CLASS(symmetry, basic) + +////////// +// default constructor, destructor, copy constructor assignment operator and helpers +////////// + +symmetry::symmetry() : type(none) +{ + debugmsg("symmetry default constructor", LOGLEVEL_CONSTRUCT); + tinfo_key = TINFO_symmetry; +} + +void symmetry::copy(const symmetry & other) +{ + inherited::copy(other); + type = other.type; + indices = other.indices; + children = other.children; +} + +DEFAULT_DESTROY(symmetry) + +////////// +// other constructors +////////// + +symmetry::symmetry(unsigned i) : type(none) +{ + debugmsg("symmetry constructor from unsigned", LOGLEVEL_CONSTRUCT); + indices.insert(i); + tinfo_key = TINFO_symmetry; +} + +symmetry::symmetry(symmetry_type t, const symmetry &c1, const symmetry &c2) : type(t) +{ + debugmsg("symmetry constructor from symmetry_type,symmetry &,symmetry &", LOGLEVEL_CONSTRUCT); + add(c1); add(c2); + tinfo_key = TINFO_symmetry; +} + +////////// +// archiving +////////// + +/** Construct object from archive_node. */ +symmetry::symmetry(const archive_node &n, const lst &sym_lst) : inherited(n, sym_lst) +{ + debugmsg("symmetry ctor from archive_node", LOGLEVEL_CONSTRUCT); + + unsigned t; + if (!(n.find_unsigned("type", t))) + throw (std::runtime_error("unknown symmetry type in archive")); + type = (symmetry_type)t; + + unsigned i = 0; + while (true) { + ex e; + if (n.find_ex("child", e, sym_lst, i)) + add(ex_to_symmetry(e)); + else + break; + i++; + } + + if (i == 0) { + while (true) { + unsigned u; + if (n.find_unsigned("index", u, i)) + indices.insert(u); + else + break; + i++; + } + } +} + +/** Archive the object. */ +void symmetry::archive(archive_node &n) const +{ + inherited::archive(n); + + n.add_unsigned("type", type); + + if (children.empty()) { + std::set::const_iterator i = indices.begin(), iend = indices.end(); + while (i != iend) { + n.add_unsigned("index", *i); + i++; + } + } else { + exvector::const_iterator i = children.begin(), iend = children.end(); + while (i != iend) { + n.add_ex("child", *i); + i++; + } + } +} + +DEFAULT_UNARCHIVE(symmetry) + +////////// +// functions overriding virtual functions from bases classes +////////// + +int symmetry::compare_same_type(const basic & other) const +{ + GINAC_ASSERT(is_of_type(other, symmetry)); + const symmetry &o = static_cast(other); + + // All symmetry trees are equal. They are not supposed to appear in + // ordinary expressions anyway... + return 0; +} + +void symmetry::print(const print_context & c, unsigned level = 0) const +{ + debugmsg("symmetry print", LOGLEVEL_PRINT); + + if (children.empty()) { + if (indices.size() > 0) + c.s << *(indices.begin()); + } else { + switch (type) { + case none: c.s << '!'; break; + case symmetric: c.s << '+'; break; + case antisymmetric: c.s << '-'; break; + case cyclic: c.s << '@'; break; + default: c.s << '?'; break; + } + c.s << '('; + for (unsigned i=0; i un; + set_union(indices.begin(), indices.end(), c.indices.begin(), c.indices.end(), inserter(un, un.begin())); + if (un.size() != indices.size() + c.indices.size()) + throw (std::logic_error("symmetry::add(): the same index appears in more than one child")); + + // Set new index set + indices.swap(un); + + // Add child node + children.push_back(c); + return *this; +} + +void symmetry::validate(unsigned n) +{ + if (indices.upper_bound(n - 1) != indices.end()) + throw (std::range_error("symmetry::verify(): index values are out of range")); + if (type != none && indices.empty()) { + for (unsigned i=0; i { + exvector::iterator v; + +public: + sy_is_less(exvector::iterator v_) : v(v_) {} + + bool operator() (const ex &lh, const ex &rh) const + { + GINAC_ASSERT(is_ex_exactly_of_type(lh, symmetry)); + GINAC_ASSERT(is_ex_exactly_of_type(rh, symmetry)); + GINAC_ASSERT(ex_to_symmetry(lh).indices.size() == ex_to_symmetry(rh).indices.size()); + std::set::const_iterator ait = ex_to_symmetry(lh).indices.begin(), aitend = ex_to_symmetry(lh).indices.end(), bit = ex_to_symmetry(rh).indices.begin(); + while (ait != aitend) { + int cmpval = v[*ait].compare(v[*bit]); + if (cmpval < 0) + return true; + else if (cmpval > 0) + return false; + ++ait; ++bit; + } + return false; + } +}; + +class sy_swap : public std::binary_function { + exvector::iterator v; + +public: + bool &swapped; + + sy_swap(exvector::iterator v_, bool &s) : v(v_), swapped(s) {} + + void operator() (const ex &lh, const ex &rh) + { + GINAC_ASSERT(is_ex_exactly_of_type(lh, symmetry)); + GINAC_ASSERT(is_ex_exactly_of_type(rh, symmetry)); + GINAC_ASSERT(ex_to_symmetry(lh).indices.size() == ex_to_symmetry(rh).indices.size()); + std::set::const_iterator ait = ex_to_symmetry(lh).indices.begin(), aitend = ex_to_symmetry(lh).indices.end(), bit = ex_to_symmetry(rh).indices.begin(); + while (ait != aitend) { + v[*ait].swap(v[*bit]); + ++ait; ++bit; + } + swapped = true; + } +}; + +int canonicalize(exvector::iterator v, const symmetry &symm) +{ + // No children? Then do nothing + if (symm.children.empty()) + return INT_MAX; + + // Canonicalize children first + bool something_changed = false; + int sign = 1; + exvector::const_iterator first = symm.children.begin(), last = symm.children.end(); + while (first != last) { + GINAC_ASSERT(is_ex_exactly_of_type(*first, symmetry)); + int child_sign = canonicalize(v, ex_to_symmetry(*first)); + if (child_sign == 0) + return 0; + if (child_sign != INT_MAX) { + something_changed = true; + sign *= child_sign; + } + first++; + } + + // Now reorder the children + first = symm.children.begin(); + switch (symm.type) { + case symmetry::symmetric: + shaker_sort(first, last, sy_is_less(v), sy_swap(v, something_changed)); + break; + case symmetry::antisymmetric: + sign *= permutation_sign(first, last, sy_is_less(v), sy_swap(v, something_changed)); + break; + case symmetry::cyclic: + cyclic_permutation(first, last, min_element(first, last, sy_is_less(v)), sy_swap(v, something_changed)); + break; + default: + break; + } + return something_changed ? sign : INT_MAX; +} + + +// Symmetrize/antisymmetrize over a vector of objects +static ex symm(const ex & e, exvector::const_iterator first, exvector::const_iterator last, bool asymmetric) +{ + // Need at least 2 objects for this operation + int num = last - first; + if (num < 2) + return e; + + // Transform object vector to a list + exlist iv_lst; + iv_lst.insert(iv_lst.begin(), first, last); + lst orig_lst(iv_lst, true); + + // Create index vectors for permutation + unsigned *iv = new unsigned[num], *iv2; + for (unsigned i=0; i + +#include "ex.h" + +namespace GiNaC { + + +class sy_is_less; +class sy_swap; + +/** This class describes the symmetry of a group of indices. These objects + * can be grouped into a tree to form complex mixed symmetries. */ +class symmetry : public basic +{ + friend class sy_is_less; + friend class sy_swap; + friend int canonicalize(exvector::iterator v, const symmetry &symm); + + GINAC_DECLARE_REGISTERED_CLASS(symmetry, basic) + + // types +public: + /** Type of symmetry */ + typedef enum { + none, /**< no symmetry properties */ + symmetric, /**< totally symmetric */ + antisymmetric, /**< totally antisymmetric */ + cyclic /**< cyclic symmetry */ + } symmetry_type; + + // other constructors +public: + /** Create leaf node that represents one index. */ + symmetry(unsigned i); + + /** Create node with two children. */ + symmetry(symmetry_type t, const symmetry &c1, const symmetry &c2); + + // functions overriding virtual functions from base classes +public: + void print(const print_context & c, unsigned level = 0) const; + + // non-virtual functions in this class +public: + /** Get symmetry type. */ + symmetry_type get_type() const {return type;} + + /** Set symmetry type. */ + void set_type(symmetry_type t) {type = t;} + + /** Add child node, check index sets for consistency. */ + symmetry &add(const symmetry &c); + + /** Verify that all indices of this node are in the range [0..n-1]. + * This function throws an exception if the verification fails. + * If the top node has a type != none and no children, add all indices + * in the range [0..n-1] as children. */ + void validate(unsigned n); + + /** Check whether this node actually represents any kind of symmetry. */ + bool has_symmetry(void) const {return type != none || !children.empty(); } + + // member variables +private: + /** Type of symmetry described by this node. */ + symmetry_type type; + + /** Sorted union set of all indices handled by this node. */ + std::set indices; + + /** Vector of child nodes. */ + exvector children; +}; + + +// global functions +inline const symmetry &ex_to_symmetry(const ex &e) +{ + return static_cast(*e.bp); +} + +inline symmetry &ex_to_nonconst_symmetry(const ex &e) +{ + return static_cast(*e.bp); +} + +inline symmetry sy_none(void) { return symmetry(); } +inline symmetry sy_none(const symmetry &c1, const symmetry &c2) { return symmetry(symmetry::none, c1, c2); } +inline symmetry sy_none(const symmetry &c1, const symmetry &c2, const symmetry &c3) { return symmetry(symmetry::none, c1, c2).add(c3); } +inline symmetry sy_none(const symmetry &c1, const symmetry &c2, const symmetry &c3, const symmetry &c4) { return symmetry(symmetry::none, c1, c2).add(c3).add(c4); } + +inline symmetry sy_symm(void) { symmetry s; s.set_type(symmetry::symmetric); return s; } +inline symmetry sy_symm(const symmetry &c1, const symmetry &c2) { return symmetry(symmetry::symmetric, c1, c2); } +inline symmetry sy_symm(const symmetry &c1, const symmetry &c2, const symmetry &c3) { return symmetry(symmetry::symmetric, c1, c2).add(c3); } +inline symmetry sy_symm(const symmetry &c1, const symmetry &c2, const symmetry &c3, const symmetry &c4) { return symmetry(symmetry::symmetric, c1, c2).add(c3).add(c4); } + +inline symmetry sy_anti(void) { symmetry s; s.set_type(symmetry::antisymmetric); return s; } +inline symmetry sy_anti(const symmetry &c1, const symmetry &c2) { return symmetry(symmetry::antisymmetric, c1, c2); } +inline symmetry sy_anti(const symmetry &c1, const symmetry &c2, const symmetry &c3) { return symmetry(symmetry::antisymmetric, c1, c2).add(c3); } +inline symmetry sy_anti(const symmetry &c1, const symmetry &c2, const symmetry &c3, const symmetry &c4) { return symmetry(symmetry::antisymmetric, c1, c2).add(c3).add(c4); } + +inline symmetry sy_cycl(void) { symmetry s; s.set_type(symmetry::cyclic); return s; } +inline symmetry sy_cycl(const symmetry &c1, const symmetry &c2) { return symmetry(symmetry::cyclic, c1, c2); } +inline symmetry sy_cycl(const symmetry &c1, const symmetry &c2, const symmetry &c3) { return symmetry(symmetry::cyclic, c1, c2).add(c3); } +inline symmetry sy_cycl(const symmetry &c1, const symmetry &c2, const symmetry &c3, const symmetry &c4) { return symmetry(symmetry::cyclic, c1, c2).add(c3).add(c4); } + +/** Canonicalize the order of elements of an expression vector, according to + * the symmetry properties defined in a symmetry tree. + * + * @param v Start of expression vector + * @param symm Root node of symmetry tree + * @return the overall sign introduced by the reordering (+1, -1 or 0) + * or INT_MAX if nothing changed */ +extern int canonicalize(exvector::iterator v, const symmetry &symm); + +/** Symmetrize expression over a set of objects (symbols, indices). */ +ex symmetrize(const ex & e, exvector::const_iterator first, exvector::const_iterator last); + +/** Symmetrize expression over a set of objects (symbols, indices). */ +inline ex symmetrize(const ex & e, const exvector & v) +{ + return symmetrize(e, v.begin(), v.end()); +} + +/** Antisymmetrize expression over a set of objects (symbols, indices). */ +ex antisymmetrize(const ex & e, exvector::const_iterator first, exvector::const_iterator last); + +/** Antisymmetrize expression over a set of objects (symbols, indices). */ +inline ex antisymmetrize(const ex & e, const exvector & v) +{ + return antisymmetrize(e, v.begin(), v.end()); +} + +/** Symmetrize expression by cyclic permuation over a set of objects + * (symbols, indices). */ +ex symmetrize_cyclic(const ex & e, exvector::const_iterator first, exvector::const_iterator last); + +/** Symmetrize expression by cyclic permutation over a set of objects + * (symbols, indices). */ +inline ex symmetrize_cyclic(const ex & e, const exvector & v) +{ + return symmetrize(e, v.begin(), v.end()); +} + +} // namespace GiNaC + +#endif // ndef __GINAC_SYMMETRY_H__