3 * Implementation of GiNaC's symmetry definitions. */
6 * GiNaC Copyright (C) 1999-2001 Johannes Gutenberg University Mainz, Germany
8 * This program is free software; you can redistribute it and/or modify
9 * it under the terms of the GNU General Public License as published by
10 * the Free Software Foundation; either version 2 of the License, or
11 * (at your option) any later version.
13 * This program is distributed in the hope that it will be useful,
14 * but WITHOUT ANY WARRANTY; without even the implied warranty of
15 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 * GNU General Public License for more details.
18 * You should have received a copy of the GNU General Public License
19 * along with this program; if not, write to the Free Software
20 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
29 #include "numeric.h" // for factorial()
37 GINAC_IMPLEMENT_REGISTERED_CLASS(symmetry, basic)
40 Some notes about the structure of a symmetry tree:
41 - The leaf nodes of the tree are of type "none", have one index, and no
42 children (of course). They are constructed by the symmetry(unsigned)
44 - Leaf nodes are the only nodes that only have one index.
45 - Container nodes contain two or more children. The "indices" set member
46 is the set union of the index sets of all children, and the "children"
47 vector stores the children themselves.
48 - The index set of each child of a "symm", "anti" or "cycl" node must
49 have the same size. It follows that the children of such a node are
50 either all leaf nodes, or all container nodes with two or more indices.
54 // default constructor, destructor, copy constructor assignment operator and helpers
57 symmetry::symmetry() : type(none)
59 debugmsg("symmetry default constructor", LOGLEVEL_CONSTRUCT);
60 tinfo_key = TINFO_symmetry;
63 void symmetry::copy(const symmetry & other)
65 inherited::copy(other);
67 indices = other.indices;
68 children = other.children;
71 DEFAULT_DESTROY(symmetry)
77 symmetry::symmetry(unsigned i) : type(none)
79 debugmsg("symmetry constructor from unsigned", LOGLEVEL_CONSTRUCT);
81 tinfo_key = TINFO_symmetry;
84 symmetry::symmetry(symmetry_type t, const symmetry &c1, const symmetry &c2) : type(t)
86 debugmsg("symmetry constructor from symmetry_type,symmetry &,symmetry &", LOGLEVEL_CONSTRUCT);
88 tinfo_key = TINFO_symmetry;
95 /** Construct object from archive_node. */
96 symmetry::symmetry(const archive_node &n, const lst &sym_lst) : inherited(n, sym_lst)
98 debugmsg("symmetry ctor from archive_node", LOGLEVEL_CONSTRUCT);
101 if (!(n.find_unsigned("type", t)))
102 throw (std::runtime_error("unknown symmetry type in archive"));
103 type = (symmetry_type)t;
108 if (n.find_ex("child", e, sym_lst, i))
109 add(ex_to<symmetry>(e));
118 if (n.find_unsigned("index", u, i))
127 /** Archive the object. */
128 void symmetry::archive(archive_node &n) const
130 inherited::archive(n);
132 n.add_unsigned("type", type);
134 if (children.empty()) {
135 std::set<unsigned>::const_iterator i = indices.begin(), iend = indices.end();
137 n.add_unsigned("index", *i);
141 exvector::const_iterator i = children.begin(), iend = children.end();
143 n.add_ex("child", *i);
149 DEFAULT_UNARCHIVE(symmetry)
152 // functions overriding virtual functions from base classes
155 int symmetry::compare_same_type(const basic & other) const
157 GINAC_ASSERT(is_of_type(other, symmetry));
159 // All symmetry trees are equal. They are not supposed to appear in
160 // ordinary expressions anyway...
164 void symmetry::print(const print_context & c, unsigned level = 0) const
166 debugmsg("symmetry print", LOGLEVEL_PRINT);
168 if (is_of_type(c, print_tree)) {
170 c.s << std::string(level, ' ') << class_name()
171 << std::hex << ", hash=0x" << hashvalue << ", flags=0x" << flags << std::dec
175 case none: c.s << "none"; break;
176 case symmetric: c.s << "symm"; break;
177 case antisymmetric: c.s << "anti"; break;
178 case cyclic: c.s << "cycl"; break;
179 default: c.s << "<unknown>"; break;
182 c.s << ", indices=(";
183 if (!indices.empty()) {
184 std::set<unsigned>::const_iterator i = indices.begin(), end = indices.end();
192 unsigned delta_indent = static_cast<const print_tree &>(c).delta_indent;
193 exvector::const_iterator i = children.begin(), end = children.end();
195 i->print(c, level + delta_indent);
201 if (children.empty()) {
202 if (indices.size() > 0)
203 c.s << *(indices.begin());
208 case none: c.s << '!'; break;
209 case symmetric: c.s << '+'; break;
210 case antisymmetric: c.s << '-'; break;
211 case cyclic: c.s << '@'; break;
212 default: c.s << '?'; break;
215 unsigned num = children.size();
216 for (unsigned i=0; i<num; i++) {
217 children[i].print(c);
227 // non-virtual functions in this class
230 symmetry &symmetry::add(const symmetry &c)
232 // All children must have the same number of indices
233 if (type != none && !children.empty()) {
234 GINAC_ASSERT(is_ex_exactly_of_type(children[0], symmetry));
235 if (ex_to<symmetry>(children[0]).indices.size() != c.indices.size())
236 throw (std::logic_error("symmetry:add(): children must have same number of indices"));
239 // Compute union of indices and check whether the two sets are disjoint
240 std::set<unsigned> un;
241 set_union(indices.begin(), indices.end(), c.indices.begin(), c.indices.end(), inserter(un, un.begin()));
242 if (un.size() != indices.size() + c.indices.size())
243 throw (std::logic_error("symmetry::add(): the same index appears in more than one child"));
249 children.push_back(c);
253 void symmetry::validate(unsigned n)
255 if (indices.upper_bound(n - 1) != indices.end())
256 throw (std::range_error("symmetry::verify(): index values are out of range"));
257 if (type != none && indices.empty()) {
258 for (unsigned i=0; i<n; i++)
267 class sy_is_less : public std::binary_function<ex, ex, bool> {
268 exvector::iterator v;
271 sy_is_less(exvector::iterator v_) : v(v_) {}
273 bool operator() (const ex &lh, const ex &rh) const
275 GINAC_ASSERT(is_ex_exactly_of_type(lh, symmetry));
276 GINAC_ASSERT(is_ex_exactly_of_type(rh, symmetry));
277 GINAC_ASSERT(ex_to<symmetry>(lh).indices.size() == ex_to<symmetry>(rh).indices.size());
278 std::set<unsigned>::const_iterator ait = ex_to<symmetry>(lh).indices.begin(), aitend = ex_to<symmetry>(lh).indices.end(), bit = ex_to<symmetry>(rh).indices.begin();
279 while (ait != aitend) {
280 int cmpval = v[*ait].compare(v[*bit]);
291 class sy_swap : public std::binary_function<ex, ex, void> {
292 exvector::iterator v;
297 sy_swap(exvector::iterator v_, bool &s) : v(v_), swapped(s) {}
299 void operator() (const ex &lh, const ex &rh)
301 GINAC_ASSERT(is_ex_exactly_of_type(lh, symmetry));
302 GINAC_ASSERT(is_ex_exactly_of_type(rh, symmetry));
303 GINAC_ASSERT(ex_to<symmetry>(lh).indices.size() == ex_to<symmetry>(rh).indices.size());
304 std::set<unsigned>::const_iterator ait = ex_to<symmetry>(lh).indices.begin(), aitend = ex_to<symmetry>(lh).indices.end(), bit = ex_to<symmetry>(rh).indices.begin();
305 while (ait != aitend) {
306 v[*ait].swap(v[*bit]);
313 int canonicalize(exvector::iterator v, const symmetry &symm)
315 // Less than two indices? Then do nothing
316 if (symm.indices.size() < 2)
319 // Canonicalize children first
320 bool something_changed = false;
322 exvector::const_iterator first = symm.children.begin(), last = symm.children.end();
323 while (first != last) {
324 GINAC_ASSERT(is_ex_exactly_of_type(*first, symmetry));
325 int child_sign = canonicalize(v, ex_to<symmetry>(*first));
328 if (child_sign != INT_MAX) {
329 something_changed = true;
335 // Now reorder the children
336 first = symm.children.begin();
338 case symmetry::symmetric:
339 // Sort the children in ascending order
340 shaker_sort(first, last, sy_is_less(v), sy_swap(v, something_changed));
342 case symmetry::antisymmetric:
343 // Sort the children in ascending order, keeping track of the signum
344 sign *= permutation_sign(first, last, sy_is_less(v), sy_swap(v, something_changed));
346 case symmetry::cyclic:
347 // Permute the smallest child to the front
348 cyclic_permutation(first, last, min_element(first, last, sy_is_less(v)), sy_swap(v, something_changed));
353 return something_changed ? sign : INT_MAX;
357 // Symmetrize/antisymmetrize over a vector of objects
358 static ex symm(const ex & e, exvector::const_iterator first, exvector::const_iterator last, bool asymmetric)
360 // Need at least 2 objects for this operation
361 unsigned num = last - first;
365 // Transform object vector to a list
367 iv_lst.insert(iv_lst.begin(), first, last);
368 lst orig_lst(iv_lst, true);
370 // Create index vectors for permutation
371 unsigned *iv = new unsigned[num], *iv2;
372 for (unsigned i=0; i<num; i++)
374 iv2 = (asymmetric ? new unsigned[num] : NULL);
376 // Loop over all permutations (the first permutation, which is the
377 // identity, is unrolled)
379 while (std::next_permutation(iv, iv + num)) {
381 for (unsigned i=0; i<num; i++)
382 new_lst.append(orig_lst.op(iv[i]));
383 ex term = e.subs(orig_lst, new_lst);
385 memcpy(iv2, iv, num * sizeof(unsigned));
386 term *= permutation_sign(iv2, iv2 + num);
394 return sum / factorial(numeric(num));
397 ex symmetrize(const ex & e, exvector::const_iterator first, exvector::const_iterator last)
399 return symm(e, first, last, false);
402 ex antisymmetrize(const ex & e, exvector::const_iterator first, exvector::const_iterator last)
404 return symm(e, first, last, true);
407 ex symmetrize_cyclic(const ex & e, exvector::const_iterator first, exvector::const_iterator last)
409 // Need at least 2 objects for this operation
410 unsigned num = last - first;
414 // Transform object vector to a list
416 iv_lst.insert(iv_lst.begin(), first, last);
417 lst orig_lst(iv_lst, true);
418 lst new_lst = orig_lst;
420 // Loop over all cyclic permutations (the first permutation, which is
421 // the identity, is unrolled)
423 for (unsigned i=0; i<num-1; i++) {
424 ex perm = new_lst.op(0);
425 new_lst.remove_first().append(perm);
426 sum += e.subs(orig_lst, new_lst);
431 /** Symmetrize expression over a list of objects (symbols, indices). */
432 ex ex::symmetrize(const lst & l) const
436 for (unsigned i=0; i<l.nops(); i++)
437 v.push_back(l.op(i));
438 return symm(*this, v.begin(), v.end(), false);
441 /** Antisymmetrize expression over a list of objects (symbols, indices). */
442 ex ex::antisymmetrize(const lst & l) const
446 for (unsigned i=0; i<l.nops(); i++)
447 v.push_back(l.op(i));
448 return symm(*this, v.begin(), v.end(), true);
451 /** Symmetrize expression by cyclic permutation over a list of objects
452 * (symbols, indices). */
453 ex ex::symmetrize_cyclic(const lst & l) const
457 for (unsigned i=0; i<l.nops(); i++)
458 v.push_back(l.op(i));
459 return GiNaC::symmetrize_cyclic(*this, v.begin(), v.end());