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()
37 GINAC_IMPLEMENT_REGISTERED_CLASS(symmetry, basic)
40 // default constructor, destructor, copy constructor assignment operator and helpers
43 symmetry::symmetry() : type(none)
45 debugmsg("symmetry default constructor", LOGLEVEL_CONSTRUCT);
46 tinfo_key = TINFO_symmetry;
49 void symmetry::copy(const symmetry & other)
51 inherited::copy(other);
53 indices = other.indices;
54 children = other.children;
57 DEFAULT_DESTROY(symmetry)
63 symmetry::symmetry(unsigned i) : type(none)
65 debugmsg("symmetry constructor from unsigned", LOGLEVEL_CONSTRUCT);
67 tinfo_key = TINFO_symmetry;
70 symmetry::symmetry(symmetry_type t, const symmetry &c1, const symmetry &c2) : type(t)
72 debugmsg("symmetry constructor from symmetry_type,symmetry &,symmetry &", LOGLEVEL_CONSTRUCT);
74 tinfo_key = TINFO_symmetry;
81 /** Construct object from archive_node. */
82 symmetry::symmetry(const archive_node &n, const lst &sym_lst) : inherited(n, sym_lst)
84 debugmsg("symmetry ctor from archive_node", LOGLEVEL_CONSTRUCT);
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 bases classes
141 int symmetry::compare_same_type(const basic & other) const
143 GINAC_ASSERT(is_of_type(other, symmetry));
144 const symmetry &o = static_cast<const symmetry &>(other);
146 // All symmetry trees are equal. They are not supposed to appear in
147 // ordinary expressions anyway...
151 void symmetry::print(const print_context & c, unsigned level = 0) const
153 debugmsg("symmetry print", LOGLEVEL_PRINT);
155 if (children.empty()) {
156 if (indices.size() > 0)
157 c.s << *(indices.begin());
160 case none: c.s << '!'; break;
161 case symmetric: c.s << '+'; break;
162 case antisymmetric: c.s << '-'; break;
163 case cyclic: c.s << '@'; break;
164 default: c.s << '?'; break;
167 for (unsigned i=0; i<children.size(); i++) {
168 children[i].print(c);
169 if (i != children.size() - 1)
177 // non-virtual functions in this class
180 symmetry &symmetry::add(const symmetry &c)
182 // All children must have the same number of indices
183 if (type != none && !children.empty()) {
184 GINAC_ASSERT(is_ex_exactly_of_type(children[0], symmetry));
185 if (ex_to<symmetry>(children[0]).indices.size() != c.indices.size())
186 throw (std::logic_error("symmetry:add(): children must have same number of indices"));
189 // Compute union of indices and check whether the two sets are disjoint
190 std::set<unsigned> un;
191 set_union(indices.begin(), indices.end(), c.indices.begin(), c.indices.end(), inserter(un, un.begin()));
192 if (un.size() != indices.size() + c.indices.size())
193 throw (std::logic_error("symmetry::add(): the same index appears in more than one child"));
199 children.push_back(c);
203 void symmetry::validate(unsigned n)
205 if (indices.upper_bound(n - 1) != indices.end())
206 throw (std::range_error("symmetry::verify(): index values are out of range"));
207 if (type != none && indices.empty()) {
208 for (unsigned i=0; i<n; i++)
217 class sy_is_less : public std::binary_function<ex, ex, bool> {
218 exvector::iterator v;
221 sy_is_less(exvector::iterator v_) : v(v_) {}
223 bool operator() (const ex &lh, const ex &rh) const
225 GINAC_ASSERT(is_ex_exactly_of_type(lh, symmetry));
226 GINAC_ASSERT(is_ex_exactly_of_type(rh, symmetry));
227 GINAC_ASSERT(ex_to<symmetry>(lh).indices.size() == ex_to<symmetry>(rh).indices.size());
228 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();
229 while (ait != aitend) {
230 int cmpval = v[*ait].compare(v[*bit]);
241 class sy_swap : public std::binary_function<ex, ex, void> {
242 exvector::iterator v;
247 sy_swap(exvector::iterator v_, bool &s) : v(v_), swapped(s) {}
249 void operator() (const ex &lh, const ex &rh)
251 GINAC_ASSERT(is_ex_exactly_of_type(lh, symmetry));
252 GINAC_ASSERT(is_ex_exactly_of_type(rh, symmetry));
253 GINAC_ASSERT(ex_to<symmetry>(lh).indices.size() == ex_to<symmetry>(rh).indices.size());
254 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();
255 while (ait != aitend) {
256 v[*ait].swap(v[*bit]);
263 int canonicalize(exvector::iterator v, const symmetry &symm)
265 // No children? Then do nothing
266 if (symm.children.empty())
269 // Canonicalize children first
270 bool something_changed = false;
272 exvector::const_iterator first = symm.children.begin(), last = symm.children.end();
273 while (first != last) {
274 GINAC_ASSERT(is_ex_exactly_of_type(*first, symmetry));
275 int child_sign = canonicalize(v, ex_to<symmetry>(*first));
278 if (child_sign != INT_MAX) {
279 something_changed = true;
285 // Now reorder the children
286 first = symm.children.begin();
288 case symmetry::symmetric:
289 // Sort the children in ascending order
290 shaker_sort(first, last, sy_is_less(v), sy_swap(v, something_changed));
292 case symmetry::antisymmetric:
293 // Sort the children in ascending order, keeping track of the signum
294 sign *= permutation_sign(first, last, sy_is_less(v), sy_swap(v, something_changed));
296 case symmetry::cyclic:
297 // Permute the smallest child to the front
298 cyclic_permutation(first, last, min_element(first, last, sy_is_less(v)), sy_swap(v, something_changed));
303 return something_changed ? sign : INT_MAX;
307 // Symmetrize/antisymmetrize over a vector of objects
308 static ex symm(const ex & e, exvector::const_iterator first, exvector::const_iterator last, bool asymmetric)
310 // Need at least 2 objects for this operation
311 int num = last - first;
315 // Transform object vector to a list
317 iv_lst.insert(iv_lst.begin(), first, last);
318 lst orig_lst(iv_lst, true);
320 // Create index vectors for permutation
321 unsigned *iv = new unsigned[num], *iv2;
322 for (unsigned i=0; i<num; i++)
324 iv2 = (asymmetric ? new unsigned[num] : NULL);
326 // Loop over all permutations (the first permutation, which is the
327 // identity, is unrolled)
329 while (std::next_permutation(iv, iv + num)) {
331 for (unsigned i=0; i<num; i++)
332 new_lst.append(orig_lst.op(iv[i]));
333 ex term = e.subs(orig_lst, new_lst);
335 memcpy(iv2, iv, num * sizeof(unsigned));
336 term *= permutation_sign(iv2, iv2 + num);
344 return sum / factorial(numeric(num));
347 ex symmetrize(const ex & e, exvector::const_iterator first, exvector::const_iterator last)
349 return symm(e, first, last, false);
352 ex antisymmetrize(const ex & e, exvector::const_iterator first, exvector::const_iterator last)
354 return symm(e, first, last, true);
357 ex symmetrize_cyclic(const ex & e, exvector::const_iterator first, exvector::const_iterator last)
359 // Need at least 2 objects for this operation
360 int num = last - first;
364 // Transform object vector to a list
366 iv_lst.insert(iv_lst.begin(), first, last);
367 lst orig_lst(iv_lst, true);
368 lst new_lst = orig_lst;
370 // Loop over all cyclic permutations (the first permutation, which is
371 // the identity, is unrolled)
373 for (unsigned i=0; i<num-1; i++) {
374 ex perm = new_lst.op(0);
375 new_lst.remove_first().append(perm);
376 sum += e.subs(orig_lst, new_lst);
381 /** Symmetrize expression over a list of objects (symbols, indices). */
382 ex ex::symmetrize(const lst & l) const
386 for (unsigned i=0; i<l.nops(); i++)
387 v.push_back(l.op(i));
388 return symm(*this, v.begin(), v.end(), false);
391 /** Antisymmetrize expression over a list of objects (symbols, indices). */
392 ex ex::antisymmetrize(const lst & l) const
396 for (unsigned i=0; i<l.nops(); i++)
397 v.push_back(l.op(i));
398 return symm(*this, v.begin(), v.end(), true);
401 /** Symmetrize expression by cyclic permutation over a list of objects
402 * (symbols, indices). */
403 ex ex::symmetrize_cyclic(const lst & l) const
407 for (unsigned i=0; i<l.nops(); i++)
408 v.push_back(l.op(i));
409 return GiNaC::symmetrize_cyclic(*this, v.begin(), v.end());