3 * Implementation of GiNaC's symmetry definitions. */
6 * GiNaC Copyright (C) 1999-2008 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
30 #include "numeric.h" // for factorial()
31 #include "operators.h"
37 GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(symmetry, basic,
38 print_func<print_context>(&symmetry::do_print).
39 print_func<print_tree>(&symmetry::do_print_tree))
42 Some notes about the structure of a symmetry tree:
43 - The leaf nodes of the tree are of type "none", have one index, and no
44 children (of course). They are constructed by the symmetry(unsigned)
46 - Leaf nodes are the only nodes that only have one index.
47 - Container nodes contain two or more children. The "indices" set member
48 is the set union of the index sets of all children, and the "children"
49 vector stores the children themselves.
50 - The index set of each child of a "symm", "anti" or "cycl" node must
51 have the same size. It follows that the children of such a node are
52 either all leaf nodes, or all container nodes with two or more indices.
56 // default constructor
59 symmetry::symmetry() : type(none)
61 setflag(status_flags::evaluated | status_flags::expanded);
68 symmetry::symmetry(unsigned i) : type(none)
71 setflag(status_flags::evaluated | status_flags::expanded);
74 symmetry::symmetry(symmetry_type t, const symmetry &c1, const symmetry &c2) : type(t)
77 setflag(status_flags::evaluated | status_flags::expanded);
84 /** Construct object from archive_node. */
85 void symmetry::read_archive(const archive_node &n, lst &sym_lst)
87 inherited::read_archive(n, sym_lst);
89 if (!(n.find_unsigned("type", t)))
90 throw (std::runtime_error("unknown symmetry type in archive"));
91 type = (symmetry_type)t;
96 if (n.find_ex("child", e, sym_lst, i))
97 add(ex_to<symmetry>(e));
106 if (n.find_unsigned("index", u, i))
114 GINAC_BIND_UNARCHIVER(symmetry);
116 /** Archive the object. */
117 void symmetry::archive(archive_node &n) const
119 inherited::archive(n);
121 n.add_unsigned("type", type);
123 if (children.empty()) {
124 std::set<unsigned>::const_iterator i = indices.begin(), iend = indices.end();
126 n.add_unsigned("index", *i);
130 exvector::const_iterator i = children.begin(), iend = children.end();
132 n.add_ex("child", *i);
139 // functions overriding virtual functions from base classes
142 int symmetry::compare_same_type(const basic & other) const
144 GINAC_ASSERT(is_a<symmetry>(other));
146 // For archiving purposes we need to have an ordering of symmetries.
147 const symmetry &othersymm = ex_to<symmetry>(other);
150 if (type > othersymm.type)
152 if (type < othersymm.type)
155 // Compare the index set.
156 size_t this_size = indices.size();
157 size_t that_size = othersymm.indices.size();
158 if (this_size > that_size)
160 if (this_size < that_size)
162 typedef std::set<unsigned>::iterator set_it;
163 set_it end = indices.end();
164 for (set_it i=indices.begin(),j=othersymm.indices.begin(); i!=end; ++i,++j) {
171 // Compare the children.
172 if (children.size() > othersymm.children.size())
174 if (children.size() < othersymm.children.size())
176 for (size_t i=0; i<children.size(); ++i) {
177 int cmpval = ex_to<symmetry>(children[i])
178 .compare_same_type(ex_to<symmetry>(othersymm.children[i]));
186 unsigned symmetry::calchash() const
188 const void* this_tinfo = (const void*)typeid(*this).name();
189 unsigned v = golden_ratio_hash((p_int)this_tinfo);
193 v ^= *(indices.begin());
195 for (exvector::const_iterator i=children.begin(); i!=children.end(); ++i)
202 if (flags & status_flags::evaluated) {
203 setflag(status_flags::hash_calculated);
210 void symmetry::do_print(const print_context & c, unsigned level) const
212 if (children.empty()) {
213 if (indices.size() > 0)
214 c.s << *(indices.begin());
219 case none: c.s << '!'; break;
220 case symmetric: c.s << '+'; break;
221 case antisymmetric: c.s << '-'; break;
222 case cyclic: c.s << '@'; break;
223 default: c.s << '?'; break;
226 size_t num = children.size();
227 for (size_t i=0; i<num; i++) {
228 children[i].print(c);
236 void symmetry::do_print_tree(const print_tree & c, unsigned level) const
238 c.s << std::string(level, ' ') << class_name() << " @" << this
239 << std::hex << ", hash=0x" << hashvalue << ", flags=0x" << flags << std::dec
243 case none: c.s << "none"; break;
244 case symmetric: c.s << "symm"; break;
245 case antisymmetric: c.s << "anti"; break;
246 case cyclic: c.s << "cycl"; break;
247 default: c.s << "<unknown>"; break;
250 c.s << ", indices=(";
251 if (!indices.empty()) {
252 std::set<unsigned>::const_iterator i = indices.begin(), end = indices.end();
260 exvector::const_iterator i = children.begin(), end = children.end();
262 i->print(c, level + c.delta_indent);
268 // non-virtual functions in this class
271 bool symmetry::has_cyclic() const
276 for (exvector::const_iterator i=children.begin(); i!=children.end(); ++i)
277 if (ex_to<symmetry>(*i).has_cyclic())
283 symmetry &symmetry::add(const symmetry &c)
285 // All children must have the same number of indices
286 if (type != none && !children.empty()) {
287 GINAC_ASSERT(is_exactly_a<symmetry>(children[0]));
288 if (ex_to<symmetry>(children[0]).indices.size() != c.indices.size())
289 throw (std::logic_error("symmetry:add(): children must have same number of indices"));
292 // Compute union of indices and check whether the two sets are disjoint
293 std::set<unsigned> un;
294 set_union(indices.begin(), indices.end(), c.indices.begin(), c.indices.end(), inserter(un, un.begin()));
295 if (un.size() != indices.size() + c.indices.size())
296 throw (std::logic_error("symmetry::add(): the same index appears in more than one child"));
302 children.push_back(c);
306 void symmetry::validate(unsigned n)
308 if (indices.upper_bound(n - 1) != indices.end())
309 throw (std::range_error("symmetry::verify(): index values are out of range"));
310 if (type != none && indices.empty()) {
311 for (unsigned i=0; i<n; i++)
320 static const symmetry & index0()
322 static ex s = (new symmetry(0))->setflag(status_flags::dynallocated);
323 return ex_to<symmetry>(s);
326 static const symmetry & index1()
328 static ex s = (new symmetry(1))->setflag(status_flags::dynallocated);
329 return ex_to<symmetry>(s);
332 static const symmetry & index2()
334 static ex s = (new symmetry(2))->setflag(status_flags::dynallocated);
335 return ex_to<symmetry>(s);
338 static const symmetry & index3()
340 static ex s = (new symmetry(3))->setflag(status_flags::dynallocated);
341 return ex_to<symmetry>(s);
344 const symmetry & not_symmetric()
346 static ex s = (new symmetry)->setflag(status_flags::dynallocated);
347 return ex_to<symmetry>(s);
350 const symmetry & symmetric2()
352 static ex s = (new symmetry(symmetry::symmetric, index0(), index1()))->setflag(status_flags::dynallocated);
353 return ex_to<symmetry>(s);
356 const symmetry & symmetric3()
358 static ex s = (new symmetry(symmetry::symmetric, index0(), index1()))->add(index2()).setflag(status_flags::dynallocated);
359 return ex_to<symmetry>(s);
362 const symmetry & symmetric4()
364 static ex s = (new symmetry(symmetry::symmetric, index0(), index1()))->add(index2()).add(index3()).setflag(status_flags::dynallocated);
365 return ex_to<symmetry>(s);
368 const symmetry & antisymmetric2()
370 static ex s = (new symmetry(symmetry::antisymmetric, index0(), index1()))->setflag(status_flags::dynallocated);
371 return ex_to<symmetry>(s);
374 const symmetry & antisymmetric3()
376 static ex s = (new symmetry(symmetry::antisymmetric, index0(), index1()))->add(index2()).setflag(status_flags::dynallocated);
377 return ex_to<symmetry>(s);
380 const symmetry & antisymmetric4()
382 static ex s = (new symmetry(symmetry::antisymmetric, index0(), index1()))->add(index2()).add(index3()).setflag(status_flags::dynallocated);
383 return ex_to<symmetry>(s);
386 class sy_is_less : public std::binary_function<ex, ex, bool> {
387 exvector::iterator v;
390 sy_is_less(exvector::iterator v_) : v(v_) {}
392 bool operator() (const ex &lh, const ex &rh) const
394 GINAC_ASSERT(is_exactly_a<symmetry>(lh));
395 GINAC_ASSERT(is_exactly_a<symmetry>(rh));
396 GINAC_ASSERT(ex_to<symmetry>(lh).indices.size() == ex_to<symmetry>(rh).indices.size());
397 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();
398 while (ait != aitend) {
399 int cmpval = v[*ait].compare(v[*bit]);
410 class sy_swap : public std::binary_function<ex, ex, void> {
411 exvector::iterator v;
416 sy_swap(exvector::iterator v_, bool &s) : v(v_), swapped(s) {}
418 void operator() (const ex &lh, const ex &rh)
420 GINAC_ASSERT(is_exactly_a<symmetry>(lh));
421 GINAC_ASSERT(is_exactly_a<symmetry>(rh));
422 GINAC_ASSERT(ex_to<symmetry>(lh).indices.size() == ex_to<symmetry>(rh).indices.size());
423 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();
424 while (ait != aitend) {
425 v[*ait].swap(v[*bit]);
432 int canonicalize(exvector::iterator v, const symmetry &symm)
434 // Less than two elements? Then do nothing
435 if (symm.indices.size() < 2)
436 return std::numeric_limits<int>::max();
438 // Canonicalize children first
439 bool something_changed = false;
441 exvector::const_iterator first = symm.children.begin(), last = symm.children.end();
442 while (first != last) {
443 GINAC_ASSERT(is_exactly_a<symmetry>(*first));
444 int child_sign = canonicalize(v, ex_to<symmetry>(*first));
447 if (child_sign != std::numeric_limits<int>::max()) {
448 something_changed = true;
454 // Now reorder the children
455 first = symm.children.begin();
457 case symmetry::symmetric:
458 // Sort the children in ascending order
459 shaker_sort(first, last, sy_is_less(v), sy_swap(v, something_changed));
461 case symmetry::antisymmetric:
462 // Sort the children in ascending order, keeping track of the signum
463 sign *= permutation_sign(first, last, sy_is_less(v), sy_swap(v, something_changed));
467 case symmetry::cyclic:
468 // Permute the smallest child to the front
469 cyclic_permutation(first, last, min_element(first, last, sy_is_less(v)), sy_swap(v, something_changed));
474 return something_changed ? sign : std::numeric_limits<int>::max();
478 // Symmetrize/antisymmetrize over a vector of objects
479 static ex symm(const ex & e, exvector::const_iterator first, exvector::const_iterator last, bool asymmetric)
481 // Need at least 2 objects for this operation
482 unsigned num = last - first;
486 // Transform object vector to a lst (for subs())
487 lst orig_lst(first, last);
489 // Create index vectors for permutation
490 unsigned *iv = new unsigned[num], *iv2;
491 for (unsigned i=0; i<num; i++)
493 iv2 = (asymmetric ? new unsigned[num] : NULL);
495 // Loop over all permutations (the first permutation, which is the
496 // identity, is unrolled)
498 while (std::next_permutation(iv, iv + num)) {
500 for (unsigned i=0; i<num; i++)
501 new_lst.append(orig_lst.op(iv[i]));
502 ex term = e.subs(orig_lst, new_lst, subs_options::no_pattern|subs_options::no_index_renaming);
504 memcpy(iv2, iv, num * sizeof(unsigned));
505 term *= permutation_sign(iv2, iv2 + num);
513 return sum / factorial(numeric(num));
516 ex symmetrize(const ex & e, exvector::const_iterator first, exvector::const_iterator last)
518 return symm(e, first, last, false);
521 ex antisymmetrize(const ex & e, exvector::const_iterator first, exvector::const_iterator last)
523 return symm(e, first, last, true);
526 ex symmetrize_cyclic(const ex & e, exvector::const_iterator first, exvector::const_iterator last)
528 // Need at least 2 objects for this operation
529 unsigned num = last - first;
533 // Transform object vector to a lst (for subs())
534 lst orig_lst(first, last);
535 lst new_lst = orig_lst;
537 // Loop over all cyclic permutations (the first permutation, which is
538 // the identity, is unrolled)
540 for (unsigned i=0; i<num-1; i++) {
541 ex perm = new_lst.op(0);
542 new_lst.remove_first().append(perm);
543 sum += e.subs(orig_lst, new_lst, subs_options::no_pattern|subs_options::no_index_renaming);
548 /** Symmetrize expression over a list of objects (symbols, indices). */
549 ex ex::symmetrize(const lst & l) const
551 exvector v(l.begin(), l.end());
552 return symm(*this, v.begin(), v.end(), false);
555 /** Antisymmetrize expression over a list of objects (symbols, indices). */
556 ex ex::antisymmetrize(const lst & l) const
558 exvector v(l.begin(), l.end());
559 return symm(*this, v.begin(), v.end(), true);
562 /** Symmetrize expression by cyclic permutation over a list of objects
563 * (symbols, indices). */
564 ex ex::symmetrize_cyclic(const lst & l) const
566 exvector v(l.begin(), l.end());
567 return GiNaC::symmetrize_cyclic(*this, v.begin(), v.end());