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()
36 GINAC_IMPLEMENT_REGISTERED_CLASS(symmetry, basic)
39 Some notes about the structure of a symmetry tree:
40 - The leaf nodes of the tree are of type "none", have one index, and no
41 children (of course). They are constructed by the symmetry(unsigned)
43 - Leaf nodes are the only nodes that only have one index.
44 - Container nodes contain two or more children. The "indices" set member
45 is the set union of the index sets of all children, and the "children"
46 vector stores the children themselves.
47 - The index set of each child of a "symm", "anti" or "cycl" node must
48 have the same size. It follows that the children of such a node are
49 either all leaf nodes, or all container nodes with two or more indices.
53 // default ctor, dtor, copy ctor, assignment operator and helpers
56 symmetry::symmetry() : type(none)
58 tinfo_key = TINFO_symmetry;
61 void symmetry::copy(const symmetry & other)
63 inherited::copy(other);
65 indices = other.indices;
66 children = other.children;
69 DEFAULT_DESTROY(symmetry)
75 symmetry::symmetry(unsigned i) : type(none)
78 tinfo_key = TINFO_symmetry;
81 symmetry::symmetry(symmetry_type t, const symmetry &c1, const symmetry &c2) : type(t)
84 tinfo_key = TINFO_symmetry;
91 /** Construct object from archive_node. */
92 symmetry::symmetry(const archive_node &n, const lst &sym_lst) : inherited(n, sym_lst)
95 if (!(n.find_unsigned("type", t)))
96 throw (std::runtime_error("unknown symmetry type in archive"));
97 type = (symmetry_type)t;
102 if (n.find_ex("child", e, sym_lst, i))
103 add(ex_to<symmetry>(e));
112 if (n.find_unsigned("index", u, i))
121 /** Archive the object. */
122 void symmetry::archive(archive_node &n) const
124 inherited::archive(n);
126 n.add_unsigned("type", type);
128 if (children.empty()) {
129 std::set<unsigned>::const_iterator i = indices.begin(), iend = indices.end();
131 n.add_unsigned("index", *i);
135 exvector::const_iterator i = children.begin(), iend = children.end();
137 n.add_ex("child", *i);
143 DEFAULT_UNARCHIVE(symmetry)
146 // functions overriding virtual functions from base classes
149 int symmetry::compare_same_type(const basic & other) const
151 GINAC_ASSERT(is_a<symmetry>(other));
153 // All symmetry trees are equal. They are not supposed to appear in
154 // ordinary expressions anyway...
158 void symmetry::print(const print_context & c, unsigned level) const
160 if (is_a<print_tree>(c)) {
162 c.s << std::string(level, ' ') << class_name()
163 << std::hex << ", hash=0x" << hashvalue << ", flags=0x" << flags << std::dec
167 case none: c.s << "none"; break;
168 case symmetric: c.s << "symm"; break;
169 case antisymmetric: c.s << "anti"; break;
170 case cyclic: c.s << "cycl"; break;
171 default: c.s << "<unknown>"; break;
174 c.s << ", indices=(";
175 if (!indices.empty()) {
176 std::set<unsigned>::const_iterator i = indices.begin(), end = indices.end();
184 unsigned delta_indent = static_cast<const print_tree &>(c).delta_indent;
185 exvector::const_iterator i = children.begin(), end = children.end();
187 i->print(c, level + delta_indent);
193 if (children.empty()) {
194 if (indices.size() > 0)
195 c.s << *(indices.begin());
200 case none: c.s << '!'; break;
201 case symmetric: c.s << '+'; break;
202 case antisymmetric: c.s << '-'; break;
203 case cyclic: c.s << '@'; break;
204 default: c.s << '?'; break;
207 unsigned num = children.size();
208 for (unsigned i=0; i<num; i++) {
209 children[i].print(c);
219 // non-virtual functions in this class
222 symmetry &symmetry::add(const symmetry &c)
224 // All children must have the same number of indices
225 if (type != none && !children.empty()) {
226 GINAC_ASSERT(is_exactly_a<symmetry>(children[0]));
227 if (ex_to<symmetry>(children[0]).indices.size() != c.indices.size())
228 throw (std::logic_error("symmetry:add(): children must have same number of indices"));
231 // Compute union of indices and check whether the two sets are disjoint
232 std::set<unsigned> un;
233 set_union(indices.begin(), indices.end(), c.indices.begin(), c.indices.end(), inserter(un, un.begin()));
234 if (un.size() != indices.size() + c.indices.size())
235 throw (std::logic_error("symmetry::add(): the same index appears in more than one child"));
241 children.push_back(c);
245 void symmetry::validate(unsigned n)
247 if (indices.upper_bound(n - 1) != indices.end())
248 throw (std::range_error("symmetry::verify(): index values are out of range"));
249 if (type != none && indices.empty()) {
250 for (unsigned i=0; i<n; i++)
259 class sy_is_less : public std::binary_function<ex, ex, bool> {
260 exvector::iterator v;
263 sy_is_less(exvector::iterator v_) : v(v_) {}
265 bool operator() (const ex &lh, const ex &rh) const
267 GINAC_ASSERT(is_exactly_a<symmetry>(lh));
268 GINAC_ASSERT(is_exactly_a<symmetry>(rh));
269 GINAC_ASSERT(ex_to<symmetry>(lh).indices.size() == ex_to<symmetry>(rh).indices.size());
270 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();
271 while (ait != aitend) {
272 int cmpval = v[*ait].compare(v[*bit]);
283 class sy_swap : public std::binary_function<ex, ex, void> {
284 exvector::iterator v;
289 sy_swap(exvector::iterator v_, bool &s) : v(v_), swapped(s) {}
291 void operator() (const ex &lh, const ex &rh)
293 GINAC_ASSERT(is_exactly_a<symmetry>(lh));
294 GINAC_ASSERT(is_exactly_a<symmetry>(rh));
295 GINAC_ASSERT(ex_to<symmetry>(lh).indices.size() == ex_to<symmetry>(rh).indices.size());
296 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();
297 while (ait != aitend) {
298 v[*ait].swap(v[*bit]);
305 int canonicalize(exvector::iterator v, const symmetry &symm)
307 // Less than two indices? Then do nothing
308 if (symm.indices.size() < 2)
311 // Canonicalize children first
312 bool something_changed = false;
314 exvector::const_iterator first = symm.children.begin(), last = symm.children.end();
315 while (first != last) {
316 GINAC_ASSERT(is_exactly_a<symmetry>(*first));
317 int child_sign = canonicalize(v, ex_to<symmetry>(*first));
320 if (child_sign != INT_MAX) {
321 something_changed = true;
327 // Now reorder the children
328 first = symm.children.begin();
330 case symmetry::symmetric:
331 // Sort the children in ascending order
332 shaker_sort(first, last, sy_is_less(v), sy_swap(v, something_changed));
334 case symmetry::antisymmetric:
335 // Sort the children in ascending order, keeping track of the signum
336 sign *= permutation_sign(first, last, sy_is_less(v), sy_swap(v, something_changed));
338 case symmetry::cyclic:
339 // Permute the smallest child to the front
340 cyclic_permutation(first, last, min_element(first, last, sy_is_less(v)), sy_swap(v, something_changed));
345 return something_changed ? sign : INT_MAX;
349 // Symmetrize/antisymmetrize over a vector of objects
350 static ex symm(const ex & e, exvector::const_iterator first, exvector::const_iterator last, bool asymmetric)
352 // Need at least 2 objects for this operation
353 unsigned num = last - first;
357 // Transform object vector to a list
359 iv_lst.insert(iv_lst.begin(), first, last);
360 lst orig_lst(iv_lst, true);
362 // Create index vectors for permutation
363 unsigned *iv = new unsigned[num], *iv2;
364 for (unsigned i=0; i<num; i++)
366 iv2 = (asymmetric ? new unsigned[num] : NULL);
368 // Loop over all permutations (the first permutation, which is the
369 // identity, is unrolled)
371 while (std::next_permutation(iv, iv + num)) {
373 for (unsigned i=0; i<num; i++)
374 new_lst.append(orig_lst.op(iv[i]));
375 ex term = e.subs(orig_lst, new_lst);
377 memcpy(iv2, iv, num * sizeof(unsigned));
378 term *= permutation_sign(iv2, iv2 + num);
386 return sum / factorial(numeric(num));
389 ex symmetrize(const ex & e, exvector::const_iterator first, exvector::const_iterator last)
391 return symm(e, first, last, false);
394 ex antisymmetrize(const ex & e, exvector::const_iterator first, exvector::const_iterator last)
396 return symm(e, first, last, true);
399 ex symmetrize_cyclic(const ex & e, exvector::const_iterator first, exvector::const_iterator last)
401 // Need at least 2 objects for this operation
402 unsigned num = last - first;
406 // Transform object vector to a list
408 iv_lst.insert(iv_lst.begin(), first, last);
409 lst orig_lst(iv_lst, true);
410 lst new_lst = orig_lst;
412 // Loop over all cyclic permutations (the first permutation, which is
413 // the identity, is unrolled)
415 for (unsigned i=0; i<num-1; i++) {
416 ex perm = new_lst.op(0);
417 new_lst.remove_first().append(perm);
418 sum += e.subs(orig_lst, new_lst);
423 /** Symmetrize expression over a list of objects (symbols, indices). */
424 ex ex::symmetrize(const lst & l) const
428 for (unsigned i=0; i<l.nops(); i++)
429 v.push_back(l.op(i));
430 return symm(*this, v.begin(), v.end(), false);
433 /** Antisymmetrize expression over a list of objects (symbols, indices). */
434 ex ex::antisymmetrize(const lst & l) const
438 for (unsigned i=0; i<l.nops(); i++)
439 v.push_back(l.op(i));
440 return symm(*this, v.begin(), v.end(), true);
443 /** Symmetrize expression by cyclic permutation over a list of objects
444 * (symbols, indices). */
445 ex ex::symmetrize_cyclic(const lst & l) const
449 for (unsigned i=0; i<l.nops(); i++)
450 v.push_back(l.op(i));
451 return GiNaC::symmetrize_cyclic(*this, v.begin(), v.end());