3 * Implementation of GiNaC's symmetry definitions. */
6 * GiNaC Copyright (C) 1999-2009 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., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
25 #include "numeric.h" // for factorial()
26 #include "operators.h"
29 #include "hash_seed.h"
38 GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(symmetry, basic,
39 print_func<print_context>(&symmetry::do_print).
40 print_func<print_tree>(&symmetry::do_print_tree))
43 Some notes about the structure of a symmetry tree:
44 - The leaf nodes of the tree are of type "none", have one index, and no
45 children (of course). They are constructed by the symmetry(unsigned)
47 - Leaf nodes are the only nodes that only have one index.
48 - Container nodes contain two or more children. The "indices" set member
49 is the set union of the index sets of all children, and the "children"
50 vector stores the children themselves.
51 - The index set of each child of a "symm", "anti" or "cycl" node must
52 have the same size. It follows that the children of such a node are
53 either all leaf nodes, or all container nodes with two or more indices.
57 // default constructor
60 symmetry::symmetry() : type(none)
62 setflag(status_flags::evaluated | status_flags::expanded);
69 symmetry::symmetry(unsigned i) : type(none)
72 setflag(status_flags::evaluated | status_flags::expanded);
75 symmetry::symmetry(symmetry_type t, const symmetry &c1, const symmetry &c2) : type(t)
78 setflag(status_flags::evaluated | status_flags::expanded);
85 /** Construct object from archive_node. */
86 void symmetry::read_archive(const archive_node &n, lst &sym_lst)
88 inherited::read_archive(n, sym_lst);
90 if (!(n.find_unsigned("type", t)))
91 throw (std::runtime_error("unknown symmetry type in archive"));
92 type = (symmetry_type)t;
97 if (n.find_ex("child", e, sym_lst, i))
98 add(ex_to<symmetry>(e));
107 if (n.find_unsigned("index", u, i))
115 GINAC_BIND_UNARCHIVER(symmetry);
117 /** Archive the object. */
118 void symmetry::archive(archive_node &n) const
120 inherited::archive(n);
122 n.add_unsigned("type", type);
124 if (children.empty()) {
125 std::set<unsigned>::const_iterator i = indices.begin(), iend = indices.end();
127 n.add_unsigned("index", *i);
131 exvector::const_iterator i = children.begin(), iend = children.end();
133 n.add_ex("child", *i);
140 // functions overriding virtual functions from base classes
143 int symmetry::compare_same_type(const basic & other) const
145 GINAC_ASSERT(is_a<symmetry>(other));
147 // For archiving purposes we need to have an ordering of symmetries.
148 const symmetry &othersymm = ex_to<symmetry>(other);
151 if (type > othersymm.type)
153 if (type < othersymm.type)
156 // Compare the index set.
157 size_t this_size = indices.size();
158 size_t that_size = othersymm.indices.size();
159 if (this_size > that_size)
161 if (this_size < that_size)
163 typedef std::set<unsigned>::iterator set_it;
164 set_it end = indices.end();
165 for (set_it i=indices.begin(),j=othersymm.indices.begin(); i!=end; ++i,++j) {
172 // Compare the children.
173 if (children.size() > othersymm.children.size())
175 if (children.size() < othersymm.children.size())
177 for (size_t i=0; i<children.size(); ++i) {
178 int cmpval = ex_to<symmetry>(children[i])
179 .compare_same_type(ex_to<symmetry>(othersymm.children[i]));
187 unsigned symmetry::calchash() const
189 unsigned v = make_hash_seed(typeid(*this));
193 if (!indices.empty())
194 v ^= *(indices.begin());
196 for (exvector::const_iterator i=children.begin(); i!=children.end(); ++i)
203 if (flags & status_flags::evaluated) {
204 setflag(status_flags::hash_calculated);
211 void symmetry::do_print(const print_context & c, unsigned level) const
213 if (children.empty()) {
214 if (indices.size() > 0)
215 c.s << *(indices.begin());
220 case none: c.s << '!'; break;
221 case symmetric: c.s << '+'; break;
222 case antisymmetric: c.s << '-'; break;
223 case cyclic: c.s << '@'; break;
224 default: c.s << '?'; break;
227 size_t num = children.size();
228 for (size_t i=0; i<num; i++) {
229 children[i].print(c);
237 void symmetry::do_print_tree(const print_tree & c, unsigned level) const
239 c.s << std::string(level, ' ') << class_name() << " @" << this
240 << std::hex << ", hash=0x" << hashvalue << ", flags=0x" << flags << std::dec
244 case none: c.s << "none"; break;
245 case symmetric: c.s << "symm"; break;
246 case antisymmetric: c.s << "anti"; break;
247 case cyclic: c.s << "cycl"; break;
248 default: c.s << "<unknown>"; break;
251 c.s << ", indices=(";
252 if (!indices.empty()) {
253 std::set<unsigned>::const_iterator i = indices.begin(), end = indices.end();
261 exvector::const_iterator i = children.begin(), end = children.end();
263 i->print(c, level + c.delta_indent);
269 // non-virtual functions in this class
272 bool symmetry::has_cyclic() const
277 for (exvector::const_iterator i=children.begin(); i!=children.end(); ++i)
278 if (ex_to<symmetry>(*i).has_cyclic())
284 symmetry &symmetry::add(const symmetry &c)
286 // All children must have the same number of indices
287 if (type != none && !children.empty()) {
288 GINAC_ASSERT(is_exactly_a<symmetry>(children[0]));
289 if (ex_to<symmetry>(children[0]).indices.size() != c.indices.size())
290 throw (std::logic_error("symmetry:add(): children must have same number of indices"));
293 // Compute union of indices and check whether the two sets are disjoint
294 std::set<unsigned> un;
295 set_union(indices.begin(), indices.end(), c.indices.begin(), c.indices.end(), inserter(un, un.begin()));
296 if (un.size() != indices.size() + c.indices.size())
297 throw (std::logic_error("symmetry::add(): the same index appears in more than one child"));
303 children.push_back(c);
307 void symmetry::validate(unsigned n)
309 if (indices.upper_bound(n - 1) != indices.end())
310 throw (std::range_error("symmetry::verify(): index values are out of range"));
311 if (type != none && indices.empty()) {
312 for (unsigned i=0; i<n; i++)
321 static const symmetry & index0()
323 static ex s = (new symmetry(0))->setflag(status_flags::dynallocated);
324 return ex_to<symmetry>(s);
327 static const symmetry & index1()
329 static ex s = (new symmetry(1))->setflag(status_flags::dynallocated);
330 return ex_to<symmetry>(s);
333 static const symmetry & index2()
335 static ex s = (new symmetry(2))->setflag(status_flags::dynallocated);
336 return ex_to<symmetry>(s);
339 static const symmetry & index3()
341 static ex s = (new symmetry(3))->setflag(status_flags::dynallocated);
342 return ex_to<symmetry>(s);
345 const symmetry & not_symmetric()
347 static ex s = (new symmetry)->setflag(status_flags::dynallocated);
348 return ex_to<symmetry>(s);
351 const symmetry & symmetric2()
353 static ex s = (new symmetry(symmetry::symmetric, index0(), index1()))->setflag(status_flags::dynallocated);
354 return ex_to<symmetry>(s);
357 const symmetry & symmetric3()
359 static ex s = (new symmetry(symmetry::symmetric, index0(), index1()))->add(index2()).setflag(status_flags::dynallocated);
360 return ex_to<symmetry>(s);
363 const symmetry & symmetric4()
365 static ex s = (new symmetry(symmetry::symmetric, index0(), index1()))->add(index2()).add(index3()).setflag(status_flags::dynallocated);
366 return ex_to<symmetry>(s);
369 const symmetry & antisymmetric2()
371 static ex s = (new symmetry(symmetry::antisymmetric, index0(), index1()))->setflag(status_flags::dynallocated);
372 return ex_to<symmetry>(s);
375 const symmetry & antisymmetric3()
377 static ex s = (new symmetry(symmetry::antisymmetric, index0(), index1()))->add(index2()).setflag(status_flags::dynallocated);
378 return ex_to<symmetry>(s);
381 const symmetry & antisymmetric4()
383 static ex s = (new symmetry(symmetry::antisymmetric, index0(), index1()))->add(index2()).add(index3()).setflag(status_flags::dynallocated);
384 return ex_to<symmetry>(s);
387 class sy_is_less : public std::binary_function<ex, ex, bool> {
388 exvector::iterator v;
391 sy_is_less(exvector::iterator v_) : v(v_) {}
393 bool operator() (const ex &lh, const ex &rh) const
395 GINAC_ASSERT(is_exactly_a<symmetry>(lh));
396 GINAC_ASSERT(is_exactly_a<symmetry>(rh));
397 GINAC_ASSERT(ex_to<symmetry>(lh).indices.size() == ex_to<symmetry>(rh).indices.size());
398 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();
399 while (ait != aitend) {
400 int cmpval = v[*ait].compare(v[*bit]);
411 class sy_swap : public std::binary_function<ex, ex, void> {
412 exvector::iterator v;
417 sy_swap(exvector::iterator v_, bool &s) : v(v_), swapped(s) {}
419 void operator() (const ex &lh, const ex &rh)
421 GINAC_ASSERT(is_exactly_a<symmetry>(lh));
422 GINAC_ASSERT(is_exactly_a<symmetry>(rh));
423 GINAC_ASSERT(ex_to<symmetry>(lh).indices.size() == ex_to<symmetry>(rh).indices.size());
424 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();
425 while (ait != aitend) {
426 v[*ait].swap(v[*bit]);
433 int canonicalize(exvector::iterator v, const symmetry &symm)
435 // Less than two elements? Then do nothing
436 if (symm.indices.size() < 2)
437 return std::numeric_limits<int>::max();
439 // Canonicalize children first
440 bool something_changed = false;
442 exvector::const_iterator first = symm.children.begin(), last = symm.children.end();
443 while (first != last) {
444 GINAC_ASSERT(is_exactly_a<symmetry>(*first));
445 int child_sign = canonicalize(v, ex_to<symmetry>(*first));
448 if (child_sign != std::numeric_limits<int>::max()) {
449 something_changed = true;
455 // Now reorder the children
456 first = symm.children.begin();
458 case symmetry::symmetric:
459 // Sort the children in ascending order
460 shaker_sort(first, last, sy_is_less(v), sy_swap(v, something_changed));
462 case symmetry::antisymmetric:
463 // Sort the children in ascending order, keeping track of the signum
464 sign *= permutation_sign(first, last, sy_is_less(v), sy_swap(v, something_changed));
468 case symmetry::cyclic:
469 // Permute the smallest child to the front
470 cyclic_permutation(first, last, min_element(first, last, sy_is_less(v)), sy_swap(v, something_changed));
475 return something_changed ? sign : std::numeric_limits<int>::max();
479 // Symmetrize/antisymmetrize over a vector of objects
480 static ex symm(const ex & e, exvector::const_iterator first, exvector::const_iterator last, bool asymmetric)
482 // Need at least 2 objects for this operation
483 unsigned num = last - first;
487 // Transform object vector to a lst (for subs())
488 lst orig_lst(first, last);
490 // Create index vectors for permutation
491 unsigned *iv = new unsigned[num], *iv2;
492 for (unsigned i=0; i<num; i++)
494 iv2 = (asymmetric ? new unsigned[num] : NULL);
496 // Loop over all permutations (the first permutation, which is the
497 // identity, is unrolled)
499 while (std::next_permutation(iv, iv + num)) {
501 for (unsigned i=0; i<num; i++)
502 new_lst.append(orig_lst.op(iv[i]));
503 ex term = e.subs(orig_lst, new_lst, subs_options::no_pattern|subs_options::no_index_renaming);
505 memcpy(iv2, iv, num * sizeof(unsigned));
506 term *= permutation_sign(iv2, iv2 + num);
514 return sum / factorial(numeric(num));
517 ex symmetrize(const ex & e, exvector::const_iterator first, exvector::const_iterator last)
519 return symm(e, first, last, false);
522 ex antisymmetrize(const ex & e, exvector::const_iterator first, exvector::const_iterator last)
524 return symm(e, first, last, true);
527 ex symmetrize_cyclic(const ex & e, exvector::const_iterator first, exvector::const_iterator last)
529 // Need at least 2 objects for this operation
530 unsigned num = last - first;
534 // Transform object vector to a lst (for subs())
535 lst orig_lst(first, last);
536 lst new_lst = orig_lst;
538 // Loop over all cyclic permutations (the first permutation, which is
539 // the identity, is unrolled)
541 for (unsigned i=0; i<num-1; i++) {
542 ex perm = new_lst.op(0);
543 new_lst.remove_first().append(perm);
544 sum += e.subs(orig_lst, new_lst, subs_options::no_pattern|subs_options::no_index_renaming);
549 /** Symmetrize expression over a list of objects (symbols, indices). */
550 ex ex::symmetrize(const lst & l) const
552 exvector v(l.begin(), l.end());
553 return symm(*this, v.begin(), v.end(), false);
556 /** Antisymmetrize expression over a list of objects (symbols, indices). */
557 ex ex::antisymmetrize(const lst & l) const
559 exvector v(l.begin(), l.end());
560 return symm(*this, v.begin(), v.end(), true);
563 /** Symmetrize expression by cyclic permutation over a list of objects
564 * (symbols, indices). */
565 ex ex::symmetrize_cyclic(const lst & l) const
567 exvector v(l.begin(), l.end());
568 return GiNaC::symmetrize_cyclic(*this, v.begin(), v.end());