3 * Implementation of GiNaC's symmetry definitions. */
6 * GiNaC Copyright (C) 1999-2003 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()
30 #include "operators.h"
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
57 symmetry::symmetry() : type(none)
59 tinfo_key = TINFO_symmetry;
66 symmetry::symmetry(unsigned i) : type(none)
69 tinfo_key = TINFO_symmetry;
72 symmetry::symmetry(symmetry_type t, const symmetry &c1, const symmetry &c2) : type(t)
75 tinfo_key = TINFO_symmetry;
82 /** Construct object from archive_node. */
83 symmetry::symmetry(const archive_node &n, lst &sym_lst) : inherited(n, sym_lst)
86 if (!(n.find_unsigned("type", t)))
87 throw (std::runtime_error("unknown symmetry type in archive"));
88 type = (symmetry_type)t;
93 if (n.find_ex("child", e, sym_lst, i))
94 add(ex_to<symmetry>(e));
103 if (n.find_unsigned("index", u, i))
112 /** Archive the object. */
113 void symmetry::archive(archive_node &n) const
115 inherited::archive(n);
117 n.add_unsigned("type", type);
119 if (children.empty()) {
120 std::set<unsigned>::const_iterator i = indices.begin(), iend = indices.end();
122 n.add_unsigned("index", *i);
126 exvector::const_iterator i = children.begin(), iend = children.end();
128 n.add_ex("child", *i);
134 DEFAULT_UNARCHIVE(symmetry)
137 // functions overriding virtual functions from base classes
140 int symmetry::compare_same_type(const basic & other) const
142 GINAC_ASSERT(is_a<symmetry>(other));
144 // All symmetry trees are equal. They are not supposed to appear in
145 // ordinary expressions anyway...
149 void symmetry::print(const print_context & c, unsigned level) const
151 if (is_a<print_tree>(c)) {
153 c.s << std::string(level, ' ') << class_name()
154 << std::hex << ", hash=0x" << hashvalue << ", flags=0x" << flags << std::dec
158 case none: c.s << "none"; break;
159 case symmetric: c.s << "symm"; break;
160 case antisymmetric: c.s << "anti"; break;
161 case cyclic: c.s << "cycl"; break;
162 default: c.s << "<unknown>"; break;
165 c.s << ", indices=(";
166 if (!indices.empty()) {
167 std::set<unsigned>::const_iterator i = indices.begin(), end = indices.end();
175 unsigned delta_indent = static_cast<const print_tree &>(c).delta_indent;
176 exvector::const_iterator i = children.begin(), end = children.end();
178 i->print(c, level + delta_indent);
184 if (children.empty()) {
185 if (indices.size() > 0)
186 c.s << *(indices.begin());
191 case none: c.s << '!'; break;
192 case symmetric: c.s << '+'; break;
193 case antisymmetric: c.s << '-'; break;
194 case cyclic: c.s << '@'; break;
195 default: c.s << '?'; break;
198 size_t num = children.size();
199 for (size_t i=0; i<num; i++) {
200 children[i].print(c);
210 // non-virtual functions in this class
213 symmetry &symmetry::add(const symmetry &c)
215 // All children must have the same number of indices
216 if (type != none && !children.empty()) {
217 GINAC_ASSERT(is_exactly_a<symmetry>(children[0]));
218 if (ex_to<symmetry>(children[0]).indices.size() != c.indices.size())
219 throw (std::logic_error("symmetry:add(): children must have same number of indices"));
222 // Compute union of indices and check whether the two sets are disjoint
223 std::set<unsigned> un;
224 set_union(indices.begin(), indices.end(), c.indices.begin(), c.indices.end(), inserter(un, un.begin()));
225 if (un.size() != indices.size() + c.indices.size())
226 throw (std::logic_error("symmetry::add(): the same index appears in more than one child"));
232 children.push_back(c);
236 void symmetry::validate(unsigned n)
238 if (indices.upper_bound(n - 1) != indices.end())
239 throw (std::range_error("symmetry::verify(): index values are out of range"));
240 if (type != none && indices.empty()) {
241 for (unsigned i=0; i<n; i++)
250 class sy_is_less : public std::binary_function<ex, ex, bool> {
251 exvector::iterator v;
254 sy_is_less(exvector::iterator v_) : v(v_) {}
256 bool operator() (const ex &lh, const ex &rh) const
258 GINAC_ASSERT(is_exactly_a<symmetry>(lh));
259 GINAC_ASSERT(is_exactly_a<symmetry>(rh));
260 GINAC_ASSERT(ex_to<symmetry>(lh).indices.size() == ex_to<symmetry>(rh).indices.size());
261 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();
262 while (ait != aitend) {
263 int cmpval = v[*ait].compare(v[*bit]);
274 class sy_swap : public std::binary_function<ex, ex, void> {
275 exvector::iterator v;
280 sy_swap(exvector::iterator v_, bool &s) : v(v_), swapped(s) {}
282 void operator() (const ex &lh, const ex &rh)
284 GINAC_ASSERT(is_exactly_a<symmetry>(lh));
285 GINAC_ASSERT(is_exactly_a<symmetry>(rh));
286 GINAC_ASSERT(ex_to<symmetry>(lh).indices.size() == ex_to<symmetry>(rh).indices.size());
287 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();
288 while (ait != aitend) {
289 v[*ait].swap(v[*bit]);
296 int canonicalize(exvector::iterator v, const symmetry &symm)
298 // Less than two elements? Then do nothing
299 if (symm.indices.size() < 2)
302 // Canonicalize children first
303 bool something_changed = false;
305 exvector::const_iterator first = symm.children.begin(), last = symm.children.end();
306 while (first != last) {
307 GINAC_ASSERT(is_exactly_a<symmetry>(*first));
308 int child_sign = canonicalize(v, ex_to<symmetry>(*first));
311 if (child_sign != INT_MAX) {
312 something_changed = true;
318 // Now reorder the children
319 first = symm.children.begin();
321 case symmetry::symmetric:
322 // Sort the children in ascending order
323 shaker_sort(first, last, sy_is_less(v), sy_swap(v, something_changed));
325 case symmetry::antisymmetric:
326 // Sort the children in ascending order, keeping track of the signum
327 sign *= permutation_sign(first, last, sy_is_less(v), sy_swap(v, something_changed));
331 case symmetry::cyclic:
332 // Permute the smallest child to the front
333 cyclic_permutation(first, last, min_element(first, last, sy_is_less(v)), sy_swap(v, something_changed));
338 return something_changed ? sign : INT_MAX;
342 // Symmetrize/antisymmetrize over a vector of objects
343 static ex symm(const ex & e, exvector::const_iterator first, exvector::const_iterator last, bool asymmetric)
345 // Need at least 2 objects for this operation
346 unsigned num = last - first;
350 // Transform object vector to a lst (for subs())
351 lst orig_lst(first, last);
353 // Create index vectors for permutation
354 unsigned *iv = new unsigned[num], *iv2;
355 for (unsigned i=0; i<num; i++)
357 iv2 = (asymmetric ? new unsigned[num] : NULL);
359 // Loop over all permutations (the first permutation, which is the
360 // identity, is unrolled)
362 while (std::next_permutation(iv, iv + num)) {
364 for (unsigned i=0; i<num; i++)
365 new_lst.append(orig_lst.op(iv[i]));
366 ex term = e.subs(orig_lst, new_lst, subs_options::no_pattern);
368 memcpy(iv2, iv, num * sizeof(unsigned));
369 term *= permutation_sign(iv2, iv2 + num);
377 return sum / factorial(numeric(num));
380 ex symmetrize(const ex & e, exvector::const_iterator first, exvector::const_iterator last)
382 return symm(e, first, last, false);
385 ex antisymmetrize(const ex & e, exvector::const_iterator first, exvector::const_iterator last)
387 return symm(e, first, last, true);
390 ex symmetrize_cyclic(const ex & e, exvector::const_iterator first, exvector::const_iterator last)
392 // Need at least 2 objects for this operation
393 unsigned num = last - first;
397 // Transform object vector to a lst (for subs())
398 lst orig_lst(first, last);
399 lst new_lst = orig_lst;
401 // Loop over all cyclic permutations (the first permutation, which is
402 // the identity, is unrolled)
404 for (unsigned i=0; i<num-1; i++) {
405 ex perm = new_lst.op(0);
406 new_lst.remove_first().append(perm);
407 sum += e.subs(orig_lst, new_lst, subs_options::no_pattern);
412 /** Symmetrize expression over a list of objects (symbols, indices). */
413 ex ex::symmetrize(const lst & l) const
415 exvector v(l.begin(), l.end());
416 return symm(*this, v.begin(), v.end(), false);
419 /** Antisymmetrize expression over a list of objects (symbols, indices). */
420 ex ex::antisymmetrize(const lst & l) const
422 exvector v(l.begin(), l.end());
423 return symm(*this, v.begin(), v.end(), true);
426 /** Symmetrize expression by cyclic permutation over a list of objects
427 * (symbols, indices). */
428 ex ex::symmetrize_cyclic(const lst & l) const
430 exvector v(l.begin(), l.end());
431 return GiNaC::symmetrize_cyclic(*this, v.begin(), v.end());