3 * Implementation of GiNaC's symmetry definitions. */
6 * GiNaC Copyright (C) 1999-2005 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"
36 GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(symmetry, basic,
37 print_func<print_context>(&symmetry::do_print).
38 print_func<print_tree>(&symmetry::do_print_tree))
41 Some notes about the structure of a symmetry tree:
42 - The leaf nodes of the tree are of type "none", have one index, and no
43 children (of course). They are constructed by the symmetry(unsigned)
45 - Leaf nodes are the only nodes that only have one index.
46 - Container nodes contain two or more children. The "indices" set member
47 is the set union of the index sets of all children, and the "children"
48 vector stores the children themselves.
49 - The index set of each child of a "symm", "anti" or "cycl" node must
50 have the same size. It follows that the children of such a node are
51 either all leaf nodes, or all container nodes with two or more indices.
55 // default constructor
58 symmetry::symmetry() : inherited(TINFO_symmetry), type(none)
60 setflag(status_flags::evaluated | status_flags::expanded);
67 symmetry::symmetry(unsigned i) : inherited(TINFO_symmetry), type(none)
70 setflag(status_flags::evaluated | status_flags::expanded);
73 symmetry::symmetry(symmetry_type t, const symmetry &c1, const symmetry &c2) : inherited(TINFO_symmetry), type(t)
76 setflag(status_flags::evaluated | status_flags::expanded);
83 /** Construct object from archive_node. */
84 symmetry::symmetry(const archive_node &n, lst &sym_lst) : inherited(n, sym_lst)
87 if (!(n.find_unsigned("type", t)))
88 throw (std::runtime_error("unknown symmetry type in archive"));
89 type = (symmetry_type)t;
94 if (n.find_ex("child", e, sym_lst, i))
95 add(ex_to<symmetry>(e));
104 if (n.find_unsigned("index", u, i))
113 /** Archive the object. */
114 void symmetry::archive(archive_node &n) const
116 inherited::archive(n);
118 n.add_unsigned("type", type);
120 if (children.empty()) {
121 std::set<unsigned>::const_iterator i = indices.begin(), iend = indices.end();
123 n.add_unsigned("index", *i);
127 exvector::const_iterator i = children.begin(), iend = children.end();
129 n.add_ex("child", *i);
135 DEFAULT_UNARCHIVE(symmetry)
138 // functions overriding virtual functions from base classes
141 int symmetry::compare_same_type(const basic & other) const
143 GINAC_ASSERT(is_a<symmetry>(other));
145 // All symmetry trees are equal. They are not supposed to appear in
146 // ordinary expressions anyway...
150 void symmetry::do_print(const print_context & c, unsigned level) const
152 if (children.empty()) {
153 if (indices.size() > 0)
154 c.s << *(indices.begin());
159 case none: c.s << '!'; break;
160 case symmetric: c.s << '+'; break;
161 case antisymmetric: c.s << '-'; break;
162 case cyclic: c.s << '@'; break;
163 default: c.s << '?'; break;
166 size_t num = children.size();
167 for (size_t i=0; i<num; i++) {
168 children[i].print(c);
176 void symmetry::do_print_tree(const print_tree & c, unsigned level) const
178 c.s << std::string(level, ' ') << class_name() << " @" << this
179 << std::hex << ", hash=0x" << hashvalue << ", flags=0x" << flags << std::dec
183 case none: c.s << "none"; break;
184 case symmetric: c.s << "symm"; break;
185 case antisymmetric: c.s << "anti"; break;
186 case cyclic: c.s << "cycl"; break;
187 default: c.s << "<unknown>"; break;
190 c.s << ", indices=(";
191 if (!indices.empty()) {
192 std::set<unsigned>::const_iterator i = indices.begin(), end = indices.end();
200 exvector::const_iterator i = children.begin(), end = children.end();
202 i->print(c, level + c.delta_indent);
208 // non-virtual functions in this class
211 symmetry &symmetry::add(const symmetry &c)
213 // All children must have the same number of indices
214 if (type != none && !children.empty()) {
215 GINAC_ASSERT(is_exactly_a<symmetry>(children[0]));
216 if (ex_to<symmetry>(children[0]).indices.size() != c.indices.size())
217 throw (std::logic_error("symmetry:add(): children must have same number of indices"));
220 // Compute union of indices and check whether the two sets are disjoint
221 std::set<unsigned> un;
222 set_union(indices.begin(), indices.end(), c.indices.begin(), c.indices.end(), inserter(un, un.begin()));
223 if (un.size() != indices.size() + c.indices.size())
224 throw (std::logic_error("symmetry::add(): the same index appears in more than one child"));
230 children.push_back(c);
234 void symmetry::validate(unsigned n)
236 if (indices.upper_bound(n - 1) != indices.end())
237 throw (std::range_error("symmetry::verify(): index values are out of range"));
238 if (type != none && indices.empty()) {
239 for (unsigned i=0; i<n; i++)
248 static const symmetry & index0()
250 static ex s = (new symmetry(0))->setflag(status_flags::dynallocated);
251 return ex_to<symmetry>(s);
254 static const symmetry & index1()
256 static ex s = (new symmetry(1))->setflag(status_flags::dynallocated);
257 return ex_to<symmetry>(s);
260 static const symmetry & index2()
262 static ex s = (new symmetry(2))->setflag(status_flags::dynallocated);
263 return ex_to<symmetry>(s);
266 static const symmetry & index3()
268 static ex s = (new symmetry(3))->setflag(status_flags::dynallocated);
269 return ex_to<symmetry>(s);
272 const symmetry & not_symmetric()
274 static ex s = (new symmetry)->setflag(status_flags::dynallocated);
275 return ex_to<symmetry>(s);
278 const symmetry & symmetric2()
280 static ex s = (new symmetry(symmetry::symmetric, index0(), index1()))->setflag(status_flags::dynallocated);
281 return ex_to<symmetry>(s);
284 const symmetry & symmetric3()
286 static ex s = (new symmetry(symmetry::symmetric, index0(), index1()))->add(index2()).setflag(status_flags::dynallocated);
287 return ex_to<symmetry>(s);
290 const symmetry & symmetric4()
292 static ex s = (new symmetry(symmetry::symmetric, index0(), index1()))->add(index2()).add(index3()).setflag(status_flags::dynallocated);
293 return ex_to<symmetry>(s);
296 const symmetry & antisymmetric2()
298 static ex s = (new symmetry(symmetry::antisymmetric, index0(), index1()))->setflag(status_flags::dynallocated);
299 return ex_to<symmetry>(s);
302 const symmetry & antisymmetric3()
304 static ex s = (new symmetry(symmetry::antisymmetric, index0(), index1()))->add(index2()).setflag(status_flags::dynallocated);
305 return ex_to<symmetry>(s);
308 const symmetry & antisymmetric4()
310 static ex s = (new symmetry(symmetry::antisymmetric, index0(), index1()))->add(index2()).add(index3()).setflag(status_flags::dynallocated);
311 return ex_to<symmetry>(s);
314 class sy_is_less : public std::binary_function<ex, ex, bool> {
315 exvector::iterator v;
318 sy_is_less(exvector::iterator v_) : v(v_) {}
320 bool operator() (const ex &lh, const ex &rh) const
322 GINAC_ASSERT(is_exactly_a<symmetry>(lh));
323 GINAC_ASSERT(is_exactly_a<symmetry>(rh));
324 GINAC_ASSERT(ex_to<symmetry>(lh).indices.size() == ex_to<symmetry>(rh).indices.size());
325 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();
326 while (ait != aitend) {
327 int cmpval = v[*ait].compare(v[*bit]);
338 class sy_swap : public std::binary_function<ex, ex, void> {
339 exvector::iterator v;
344 sy_swap(exvector::iterator v_, bool &s) : v(v_), swapped(s) {}
346 void operator() (const ex &lh, const ex &rh)
348 GINAC_ASSERT(is_exactly_a<symmetry>(lh));
349 GINAC_ASSERT(is_exactly_a<symmetry>(rh));
350 GINAC_ASSERT(ex_to<symmetry>(lh).indices.size() == ex_to<symmetry>(rh).indices.size());
351 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();
352 while (ait != aitend) {
353 v[*ait].swap(v[*bit]);
360 int canonicalize(exvector::iterator v, const symmetry &symm)
362 // Less than two elements? Then do nothing
363 if (symm.indices.size() < 2)
366 // Canonicalize children first
367 bool something_changed = false;
369 exvector::const_iterator first = symm.children.begin(), last = symm.children.end();
370 while (first != last) {
371 GINAC_ASSERT(is_exactly_a<symmetry>(*first));
372 int child_sign = canonicalize(v, ex_to<symmetry>(*first));
375 if (child_sign != INT_MAX) {
376 something_changed = true;
382 // Now reorder the children
383 first = symm.children.begin();
385 case symmetry::symmetric:
386 // Sort the children in ascending order
387 shaker_sort(first, last, sy_is_less(v), sy_swap(v, something_changed));
389 case symmetry::antisymmetric:
390 // Sort the children in ascending order, keeping track of the signum
391 sign *= permutation_sign(first, last, sy_is_less(v), sy_swap(v, something_changed));
395 case symmetry::cyclic:
396 // Permute the smallest child to the front
397 cyclic_permutation(first, last, min_element(first, last, sy_is_less(v)), sy_swap(v, something_changed));
402 return something_changed ? sign : INT_MAX;
406 // Symmetrize/antisymmetrize over a vector of objects
407 static ex symm(const ex & e, exvector::const_iterator first, exvector::const_iterator last, bool asymmetric)
409 // Need at least 2 objects for this operation
410 unsigned num = last - first;
414 // Transform object vector to a lst (for subs())
415 lst orig_lst(first, last);
417 // Create index vectors for permutation
418 unsigned *iv = new unsigned[num], *iv2;
419 for (unsigned i=0; i<num; i++)
421 iv2 = (asymmetric ? new unsigned[num] : NULL);
423 // Loop over all permutations (the first permutation, which is the
424 // identity, is unrolled)
426 while (std::next_permutation(iv, iv + num)) {
428 for (unsigned i=0; i<num; i++)
429 new_lst.append(orig_lst.op(iv[i]));
430 ex term = e.subs(orig_lst, new_lst, subs_options::no_pattern);
432 memcpy(iv2, iv, num * sizeof(unsigned));
433 term *= permutation_sign(iv2, iv2 + num);
441 return sum / factorial(numeric(num));
444 ex symmetrize(const ex & e, exvector::const_iterator first, exvector::const_iterator last)
446 return symm(e, first, last, false);
449 ex antisymmetrize(const ex & e, exvector::const_iterator first, exvector::const_iterator last)
451 return symm(e, first, last, true);
454 ex symmetrize_cyclic(const ex & e, exvector::const_iterator first, exvector::const_iterator last)
456 // Need at least 2 objects for this operation
457 unsigned num = last - first;
461 // Transform object vector to a lst (for subs())
462 lst orig_lst(first, last);
463 lst new_lst = orig_lst;
465 // Loop over all cyclic permutations (the first permutation, which is
466 // the identity, is unrolled)
468 for (unsigned i=0; i<num-1; i++) {
469 ex perm = new_lst.op(0);
470 new_lst.remove_first().append(perm);
471 sum += e.subs(orig_lst, new_lst, subs_options::no_pattern);
476 /** Symmetrize expression over a list of objects (symbols, indices). */
477 ex ex::symmetrize(const lst & l) const
479 exvector v(l.begin(), l.end());
480 return symm(*this, v.begin(), v.end(), false);
483 /** Antisymmetrize expression over a list of objects (symbols, indices). */
484 ex ex::antisymmetrize(const lst & l) const
486 exvector v(l.begin(), l.end());
487 return symm(*this, v.begin(), v.end(), true);
490 /** Symmetrize expression by cyclic permutation over a list of objects
491 * (symbols, indices). */
492 ex ex::symmetrize_cyclic(const lst & l) const
494 exvector v(l.begin(), l.end());
495 return GiNaC::symmetrize_cyclic(*this, v.begin(), v.end());