3 * Implementation of GiNaC's symmetry definitions. */
6 * GiNaC Copyright (C) 1999-2011 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
26 #include "numeric.h" // for factorial()
27 #include "operators.h"
30 #include "hash_seed.h"
39 GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(symmetry, basic,
40 print_func<print_context>(&symmetry::do_print).
41 print_func<print_tree>(&symmetry::do_print_tree))
44 Some notes about the structure of a symmetry tree:
45 - The leaf nodes of the tree are of type "none", have one index, and no
46 children (of course). They are constructed by the symmetry(unsigned)
48 - Leaf nodes are the only nodes that only have one index.
49 - Container nodes contain two or more children. The "indices" set member
50 is the set union of the index sets of all children, and the "children"
51 vector stores the children themselves.
52 - The index set of each child of a "symm", "anti" or "cycl" node must
53 have the same size. It follows that the children of such a node are
54 either all leaf nodes, or all container nodes with two or more indices.
58 // default constructor
61 symmetry::symmetry() : type(none)
63 setflag(status_flags::evaluated | status_flags::expanded);
70 symmetry::symmetry(unsigned i) : type(none)
73 setflag(status_flags::evaluated | status_flags::expanded);
76 symmetry::symmetry(symmetry_type t, const symmetry &c1, const symmetry &c2) : type(t)
79 setflag(status_flags::evaluated | status_flags::expanded);
86 /** Construct object from archive_node. */
87 void symmetry::read_archive(const archive_node &n, lst &sym_lst)
89 inherited::read_archive(n, sym_lst);
91 if (!(n.find_unsigned("type", t)))
92 throw (std::runtime_error("unknown symmetry type in archive"));
93 type = (symmetry_type)t;
98 if (n.find_ex("child", e, sym_lst, i))
99 add(ex_to<symmetry>(e));
108 if (n.find_unsigned("index", u, i))
116 GINAC_BIND_UNARCHIVER(symmetry);
118 /** Archive the object. */
119 void symmetry::archive(archive_node &n) const
121 inherited::archive(n);
123 n.add_unsigned("type", type);
125 if (children.empty()) {
126 std::set<unsigned>::const_iterator i = indices.begin(), iend = indices.end();
128 n.add_unsigned("index", *i);
132 exvector::const_iterator i = children.begin(), iend = children.end();
134 n.add_ex("child", *i);
141 // functions overriding virtual functions from base classes
144 int symmetry::compare_same_type(const basic & other) const
146 GINAC_ASSERT(is_a<symmetry>(other));
148 // For archiving purposes we need to have an ordering of symmetries.
149 const symmetry &othersymm = ex_to<symmetry>(other);
152 if (type > othersymm.type)
154 if (type < othersymm.type)
157 // Compare the index set.
158 size_t this_size = indices.size();
159 size_t that_size = othersymm.indices.size();
160 if (this_size > that_size)
162 if (this_size < that_size)
164 typedef std::set<unsigned>::const_iterator set_it;
165 set_it end = indices.end();
166 for (set_it i=indices.begin(),j=othersymm.indices.begin(); i!=end; ++i,++j) {
173 // Compare the children.
174 if (children.size() > othersymm.children.size())
176 if (children.size() < othersymm.children.size())
178 for (size_t i=0; i<children.size(); ++i) {
179 int cmpval = ex_to<symmetry>(children[i])
180 .compare_same_type(ex_to<symmetry>(othersymm.children[i]));
188 unsigned symmetry::calchash() const
190 unsigned v = make_hash_seed(typeid(*this));
194 if (!indices.empty())
195 v ^= *(indices.begin());
197 for (exvector::const_iterator i=children.begin(); i!=children.end(); ++i)
204 if (flags & status_flags::evaluated) {
205 setflag(status_flags::hash_calculated);
212 void symmetry::do_print(const print_context & c, unsigned level) const
214 if (children.empty()) {
215 if (indices.size() > 0)
216 c.s << *(indices.begin());
221 case none: c.s << '!'; break;
222 case symmetric: c.s << '+'; break;
223 case antisymmetric: c.s << '-'; break;
224 case cyclic: c.s << '@'; break;
225 default: c.s << '?'; break;
228 size_t num = children.size();
229 for (size_t i=0; i<num; i++) {
230 children[i].print(c);
238 void symmetry::do_print_tree(const print_tree & c, unsigned level) const
240 c.s << std::string(level, ' ') << class_name() << " @" << this
241 << std::hex << ", hash=0x" << hashvalue << ", flags=0x" << flags << std::dec
245 case none: c.s << "none"; break;
246 case symmetric: c.s << "symm"; break;
247 case antisymmetric: c.s << "anti"; break;
248 case cyclic: c.s << "cycl"; break;
249 default: c.s << "<unknown>"; break;
252 c.s << ", indices=(";
253 if (!indices.empty()) {
254 std::set<unsigned>::const_iterator i = indices.begin(), end = indices.end();
262 exvector::const_iterator i = children.begin(), end = children.end();
264 i->print(c, level + c.delta_indent);
270 // non-virtual functions in this class
273 bool symmetry::has_cyclic() const
278 for (exvector::const_iterator i=children.begin(); i!=children.end(); ++i)
279 if (ex_to<symmetry>(*i).has_cyclic())
285 symmetry &symmetry::add(const symmetry &c)
287 // All children must have the same number of indices
288 if (type != none && !children.empty()) {
289 GINAC_ASSERT(is_exactly_a<symmetry>(children[0]));
290 if (ex_to<symmetry>(children[0]).indices.size() != c.indices.size())
291 throw (std::logic_error("symmetry:add(): children must have same number of indices"));
294 // Compute union of indices and check whether the two sets are disjoint
295 std::set<unsigned> un;
296 set_union(indices.begin(), indices.end(), c.indices.begin(), c.indices.end(), inserter(un, un.begin()));
297 if (un.size() != indices.size() + c.indices.size())
298 throw (std::logic_error("symmetry::add(): the same index appears in more than one child"));
304 children.push_back(c);
308 void symmetry::validate(unsigned n)
310 if (indices.upper_bound(n - 1) != indices.end())
311 throw (std::range_error("symmetry::verify(): index values are out of range"));
312 if (type != none && indices.empty()) {
313 for (unsigned i=0; i<n; i++)
322 static const symmetry & index0()
324 static ex s = (new symmetry(0))->setflag(status_flags::dynallocated);
325 return ex_to<symmetry>(s);
328 static const symmetry & index1()
330 static ex s = (new symmetry(1))->setflag(status_flags::dynallocated);
331 return ex_to<symmetry>(s);
334 static const symmetry & index2()
336 static ex s = (new symmetry(2))->setflag(status_flags::dynallocated);
337 return ex_to<symmetry>(s);
340 static const symmetry & index3()
342 static ex s = (new symmetry(3))->setflag(status_flags::dynallocated);
343 return ex_to<symmetry>(s);
346 const symmetry & not_symmetric()
348 static ex s = (new symmetry)->setflag(status_flags::dynallocated);
349 return ex_to<symmetry>(s);
352 const symmetry & symmetric2()
354 static ex s = (new symmetry(symmetry::symmetric, index0(), index1()))->setflag(status_flags::dynallocated);
355 return ex_to<symmetry>(s);
358 const symmetry & symmetric3()
360 static ex s = (new symmetry(symmetry::symmetric, index0(), index1()))->add(index2()).setflag(status_flags::dynallocated);
361 return ex_to<symmetry>(s);
364 const symmetry & symmetric4()
366 static ex s = (new symmetry(symmetry::symmetric, index0(), index1()))->add(index2()).add(index3()).setflag(status_flags::dynallocated);
367 return ex_to<symmetry>(s);
370 const symmetry & antisymmetric2()
372 static ex s = (new symmetry(symmetry::antisymmetric, index0(), index1()))->setflag(status_flags::dynallocated);
373 return ex_to<symmetry>(s);
376 const symmetry & antisymmetric3()
378 static ex s = (new symmetry(symmetry::antisymmetric, index0(), index1()))->add(index2()).setflag(status_flags::dynallocated);
379 return ex_to<symmetry>(s);
382 const symmetry & antisymmetric4()
384 static ex s = (new symmetry(symmetry::antisymmetric, index0(), index1()))->add(index2()).add(index3()).setflag(status_flags::dynallocated);
385 return ex_to<symmetry>(s);
388 class sy_is_less : public std::binary_function<ex, ex, bool> {
389 exvector::iterator v;
392 sy_is_less(exvector::iterator v_) : v(v_) {}
394 bool operator() (const ex &lh, const ex &rh) const
396 GINAC_ASSERT(is_exactly_a<symmetry>(lh));
397 GINAC_ASSERT(is_exactly_a<symmetry>(rh));
398 GINAC_ASSERT(ex_to<symmetry>(lh).indices.size() == ex_to<symmetry>(rh).indices.size());
399 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();
400 while (ait != aitend) {
401 int cmpval = v[*ait].compare(v[*bit]);
412 class sy_swap : public std::binary_function<ex, ex, void> {
413 exvector::iterator v;
418 sy_swap(exvector::iterator v_, bool &s) : v(v_), swapped(s) {}
420 void operator() (const ex &lh, const ex &rh)
422 GINAC_ASSERT(is_exactly_a<symmetry>(lh));
423 GINAC_ASSERT(is_exactly_a<symmetry>(rh));
424 GINAC_ASSERT(ex_to<symmetry>(lh).indices.size() == ex_to<symmetry>(rh).indices.size());
425 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();
426 while (ait != aitend) {
427 v[*ait].swap(v[*bit]);
434 int canonicalize(exvector::iterator v, const symmetry &symm)
436 // Less than two elements? Then do nothing
437 if (symm.indices.size() < 2)
438 return std::numeric_limits<int>::max();
440 // Canonicalize children first
441 bool something_changed = false;
443 exvector::const_iterator first = symm.children.begin(), last = symm.children.end();
444 while (first != last) {
445 GINAC_ASSERT(is_exactly_a<symmetry>(*first));
446 int child_sign = canonicalize(v, ex_to<symmetry>(*first));
449 if (child_sign != std::numeric_limits<int>::max()) {
450 something_changed = true;
456 // Now reorder the children
457 first = symm.children.begin();
459 case symmetry::symmetric:
460 // Sort the children in ascending order
461 shaker_sort(first, last, sy_is_less(v), sy_swap(v, something_changed));
463 case symmetry::antisymmetric:
464 // Sort the children in ascending order, keeping track of the signum
465 sign *= permutation_sign(first, last, sy_is_less(v), sy_swap(v, something_changed));
469 case symmetry::cyclic:
470 // Permute the smallest child to the front
471 cyclic_permutation(first, last, min_element(first, last, sy_is_less(v)), sy_swap(v, something_changed));
476 return something_changed ? sign : std::numeric_limits<int>::max();
480 // Symmetrize/antisymmetrize over a vector of objects
481 static ex symm(const ex & e, exvector::const_iterator first, exvector::const_iterator last, bool asymmetric)
483 // Need at least 2 objects for this operation
484 unsigned num = last - first;
488 // Transform object vector to a lst (for subs())
489 lst orig_lst(first, last);
491 // Create index vectors for permutation
492 unsigned *iv = new unsigned[num], *iv2;
493 for (unsigned i=0; i<num; i++)
495 iv2 = (asymmetric ? new unsigned[num] : NULL);
497 // Loop over all permutations (the first permutation, which is the
498 // identity, is unrolled)
501 while (std::next_permutation(iv, iv + num)) {
503 for (unsigned i=0; i<num; i++)
504 new_lst.append(orig_lst.op(iv[i]));
505 ex term = e.subs(orig_lst, new_lst, subs_options::no_pattern|subs_options::no_index_renaming);
507 memcpy(iv2, iv, num * sizeof(unsigned));
508 term *= permutation_sign(iv2, iv2 + num);
510 sum_v.push_back(term);
512 ex sum = (new add(sum_v))->setflag(status_flags::dynallocated);
517 return sum / factorial(numeric(num));
520 ex symmetrize(const ex & e, exvector::const_iterator first, exvector::const_iterator last)
522 return symm(e, first, last, false);
525 ex antisymmetrize(const ex & e, exvector::const_iterator first, exvector::const_iterator last)
527 return symm(e, first, last, true);
530 ex symmetrize_cyclic(const ex & e, exvector::const_iterator first, exvector::const_iterator last)
532 // Need at least 2 objects for this operation
533 unsigned num = last - first;
537 // Transform object vector to a lst (for subs())
538 lst orig_lst(first, last);
539 lst new_lst = orig_lst;
541 // Loop over all cyclic permutations (the first permutation, which is
542 // the identity, is unrolled)
544 for (unsigned i=0; i<num-1; i++) {
545 ex perm = new_lst.op(0);
546 new_lst.remove_first().append(perm);
547 sum += e.subs(orig_lst, new_lst, subs_options::no_pattern|subs_options::no_index_renaming);
552 /** Symmetrize expression over a list of objects (symbols, indices). */
553 ex ex::symmetrize(const lst & l) const
555 exvector v(l.begin(), l.end());
556 return symm(*this, v.begin(), v.end(), false);
559 /** Antisymmetrize expression over a list of objects (symbols, indices). */
560 ex ex::antisymmetrize(const lst & l) const
562 exvector v(l.begin(), l.end());
563 return symm(*this, v.begin(), v.end(), true);
566 /** Symmetrize expression by cyclic permutation over a list of objects
567 * (symbols, indices). */
568 ex ex::symmetrize_cyclic(const lst & l) const
570 exvector v(l.begin(), l.end());
571 return GiNaC::symmetrize_cyclic(*this, v.begin(), v.end());