ginac/symmetry.cpp

   1 /** @file symmetry.cpp
   2  *
   3  *  Implementation of GiNaC's symmetry definitions. */
   4
   5 /*
   6  *  GiNaC Copyright (C) 1999-2011 Johannes Gutenberg University Mainz, Germany
   7  *
   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.
  12  *
  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.
  17  *
  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
  21  */
  22
  23 #include "symmetry.h"
  24 #include "lst.h"
  25 #include "add.h"
  26 #include "numeric.h" // for factorial()
  27 #include "operators.h"
  28 #include "archive.h"
  29 #include "utils.h"
  30 #include "hash_seed.h"
  31
  32 #include <functional>
  33 #include <iostream>
  34 #include <limits>
  35 #include <stdexcept>
  36
  37 namespace GiNaC {
  38
  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))
  42
  43 /*
  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)
  47       constructor.
  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.
  55 */
  56
  57 //////////
  58 // default constructor
  59 //////////
  60
  61 symmetry::symmetry() :  type(none)
  62 {
  63         setflag(status_flags::evaluated | status_flags::expanded);
  64 }
  65
  66 //////////
  67 // other constructors
  68 //////////
  69
  70 symmetry::symmetry(unsigned i) :  type(none)
  71 {
  72         indices.insert(i);
  73         setflag(status_flags::evaluated | status_flags::expanded);
  74 }
  75
  76 symmetry::symmetry(symmetry_type t, const symmetry &c1, const symmetry &c2) :  type(t)
  77 {
  78         add(c1); add(c2);
  79         setflag(status_flags::evaluated | status_flags::expanded);
  80 }
  81
  82 //////////
  83 // archiving
  84 //////////
  85
  86 /** Construct object from archive_node. */
  87 void symmetry::read_archive(const archive_node &n, lst &sym_lst)
  88 {
  89         inherited::read_archive(n, sym_lst);
  90         unsigned t;
  91         if (!(n.find_unsigned("type", t)))
  92                 throw (std::runtime_error("unknown symmetry type in archive"));
  93         type = (symmetry_type)t;
  94
  95         unsigned i = 0;
  96         while (true) {
  97                 ex e;
  98                 if (n.find_ex("child", e, sym_lst, i))
  99                         add(ex_to<symmetry>(e));
 100                 else
 101                         break;
 102                 i++;
 103         }
 104
 105         if (i == 0) {
 106                 while (true) {
 107                         unsigned u;
 108                         if (n.find_unsigned("index", u, i))
 109                                 indices.insert(u);
 110                         else
 111                                 break;
 112                         i++;
 113                 }
 114         }
 115 }
 116 GINAC_BIND_UNARCHIVER(symmetry);
 117
 118 /** Archive the object. */
 119 void symmetry::archive(archive_node &n) const
 120 {
 121         inherited::archive(n);
 122
 123         n.add_unsigned("type", type);
 124
 125         if (children.empty()) {
 126                 std::set<unsigned>::const_iterator i = indices.begin(), iend = indices.end();
 127                 while (i != iend) {
 128                         n.add_unsigned("index", *i);
 129                         i++;
 130                 }
 131         } else {
 132                 exvector::const_iterator i = children.begin(), iend = children.end();
 133                 while (i != iend) {
 134                         n.add_ex("child", *i);
 135                         i++;
 136                 }
 137         }
 138 }
 139
 140 //////////
 141 // functions overriding virtual functions from base classes
 142 //////////
 143
 144 int symmetry::compare_same_type(const basic & other) const
 145 {
 146         GINAC_ASSERT(is_a<symmetry>(other));
 147
 148         // For archiving purposes we need to have an ordering of symmetries.
 149         const symmetry &othersymm = ex_to<symmetry>(other);
 150
 151         // Compare type.
 152         if (type > othersymm.type)
 153                 return 1;
 154         if (type < othersymm.type)
 155                 return -1;
 156
 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)
 161                 return 1;
 162         if (this_size < that_size)
 163                 return -1;
 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) {
 167                 if(*i < *j)
 168                         return 1;
 169                 if(*i > *j)
 170                         return -1;
 171         }
 172
 173         // Compare the children.
 174         if (children.size() > othersymm.children.size())
 175                 return 1;
 176         if (children.size() < othersymm.children.size())
 177                 return -1;
 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]));
 181                 if (cmpval)
 182                         return cmpval;
 183         }
 184
 185         return 0;
 186 }
 187
 188 unsigned symmetry::calchash() const
 189 {
 190         unsigned v = make_hash_seed(typeid(*this));
 191
 192         if (type == none) {
 193                 v = rotate_left(v);
 194                 if (!indices.empty())
 195                         v ^= *(indices.begin());
 196         } else {
 197                 for (exvector::const_iterator i=children.begin(); i!=children.end(); ++i)
 198                 {
 199                         v = rotate_left(v);
 200                         v ^= i->gethash();
 201                 }
 202         }
 203
 204         if (flags & status_flags::evaluated) {
 205                 setflag(status_flags::hash_calculated);
 206                 hashvalue = v;
 207         }
 208
 209         return v;
 210 }
 211
 212 void symmetry::do_print(const print_context & c, unsigned level) const
 213 {
 214         if (children.empty()) {
 215                 if (indices.size() > 0)
 216                         c.s << *(indices.begin());
 217                 else
 218                         c.s << "none";
 219         } else {
 220                 switch (type) {
 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;
 226                 }
 227                 c.s << '(';
 228                 size_t num = children.size();
 229                 for (size_t i=0; i<num; i++) {
 230                         children[i].print(c);
 231                         if (i != num - 1)
 232                                 c.s << ",";
 233                 }
 234                 c.s << ')';
 235         }
 236 }
 237
 238 void symmetry::do_print_tree(const print_tree & c, unsigned level) const
 239 {
 240         c.s << std::string(level, ' ') << class_name() << " @" << this
 241             << std::hex << ", hash=0x" << hashvalue << ", flags=0x" << flags << std::dec
 242             << ", type=";
 243
 244         switch (type) {
 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;
 250         }
 251
 252         c.s << ", indices=(";
 253         if (!indices.empty()) {
 254                 std::set<unsigned>::const_iterator i = indices.begin(), end = indices.end();
 255                 --end;
 256                 while (i != end)
 257                         c.s << *i++ << ",";
 258                 c.s << *i;
 259         }
 260         c.s << ")\n";
 261
 262         exvector::const_iterator i = children.begin(), end = children.end();
 263         while (i != end) {
 264                 i->print(c, level + c.delta_indent);
 265                 ++i;
 266         }
 267 }
 268
 269 //////////
 270 // non-virtual functions in this class
 271 //////////
 272
 273 bool symmetry::has_cyclic() const
 274 {
 275         if (type == cyclic)
 276                 return true;
 277
 278         for (exvector::const_iterator i=children.begin(); i!=children.end(); ++i)
 279                 if (ex_to<symmetry>(*i).has_cyclic())
 280                         return true;
 281
 282         return false;
 283 }
 284
 285 symmetry &symmetry::add(const symmetry &c)
 286 {
 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"));
 292         }
 293
 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"));
 299
 300         // Set new index set
 301         indices.swap(un);
 302
 303         // Add child node
 304         children.push_back(c);
 305         return *this;
 306 }
 307
 308 void symmetry::validate(unsigned n)
 309 {
 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++)
 314                         add(i);
 315         }
 316 }
 317
 318 //////////
 319 // global functions
 320 //////////
 321
 322 static const symmetry & index0()
 323 {
 324         static ex s = (new symmetry(0))->setflag(status_flags::dynallocated);
 325         return ex_to<symmetry>(s);
 326 }
 327
 328 static const symmetry & index1()
 329 {
 330         static ex s = (new symmetry(1))->setflag(status_flags::dynallocated);
 331         return ex_to<symmetry>(s);
 332 }
 333
 334 static const symmetry & index2()
 335 {
 336         static ex s = (new symmetry(2))->setflag(status_flags::dynallocated);
 337         return ex_to<symmetry>(s);
 338 }
 339
 340 static const symmetry & index3()
 341 {
 342         static ex s = (new symmetry(3))->setflag(status_flags::dynallocated);
 343         return ex_to<symmetry>(s);
 344 }
 345
 346 const symmetry & not_symmetric()
 347 {
 348         static ex s = (new symmetry)->setflag(status_flags::dynallocated);
 349         return ex_to<symmetry>(s);
 350 }
 351
 352 const symmetry & symmetric2()
 353 {
 354         static ex s = (new symmetry(symmetry::symmetric, index0(), index1()))->setflag(status_flags::dynallocated);
 355         return ex_to<symmetry>(s);
 356 }
 357
 358 const symmetry & symmetric3()
 359 {
 360         static ex s = (new symmetry(symmetry::symmetric, index0(), index1()))->add(index2()).setflag(status_flags::dynallocated);
 361         return ex_to<symmetry>(s);
 362 }
 363
 364 const symmetry & symmetric4()
 365 {
 366         static ex s = (new symmetry(symmetry::symmetric, index0(), index1()))->add(index2()).add(index3()).setflag(status_flags::dynallocated);
 367         return ex_to<symmetry>(s);
 368 }
 369
 370 const symmetry & antisymmetric2()
 371 {
 372         static ex s = (new symmetry(symmetry::antisymmetric, index0(), index1()))->setflag(status_flags::dynallocated);
 373         return ex_to<symmetry>(s);
 374 }
 375
 376 const symmetry & antisymmetric3()
 377 {
 378         static ex s = (new symmetry(symmetry::antisymmetric, index0(), index1()))->add(index2()).setflag(status_flags::dynallocated);
 379         return ex_to<symmetry>(s);
 380 }
 381
 382 const symmetry & antisymmetric4()
 383 {
 384         static ex s = (new symmetry(symmetry::antisymmetric, index0(), index1()))->add(index2()).add(index3()).setflag(status_flags::dynallocated);
 385         return ex_to<symmetry>(s);
 386 }
 387
 388 class sy_is_less : public std::binary_function<ex, ex, bool> {
 389         exvector::iterator v;
 390
 391 public:
 392         sy_is_less(exvector::iterator v_) : v(v_) {}
 393
 394         bool operator() (const ex &lh, const ex &rh) const
 395         {
 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]);
 402                         if (cmpval < 0)
 403                                 return true;
 404                         else if (cmpval > 0)
 405                                 return false;
 406                         ++ait; ++bit;
 407                 }
 408                 return false;
 409         }
 410 };
 411
 412 class sy_swap : public std::binary_function<ex, ex, void> {
 413         exvector::iterator v;
 414
 415 public:
 416         bool &swapped;
 417
 418         sy_swap(exvector::iterator v_, bool &s) : v(v_), swapped(s) {}
 419
 420         void operator() (const ex &lh, const ex &rh)
 421         {
 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]);
 428                         ++ait; ++bit;
 429                 }
 430                 swapped = true;
 431         }
 432 };
 433
 434 int canonicalize(exvector::iterator v, const symmetry &symm)
 435 {
 436         // Less than two elements? Then do nothing
 437         if (symm.indices.size() < 2)
 438                 return std::numeric_limits<int>::max();
 439
 440         // Canonicalize children first
 441         bool something_changed = false;
 442         int sign = 1;
 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));
 447                 if (child_sign == 0)
 448                         return 0;
 449                 if (child_sign != std::numeric_limits<int>::max()) {
 450                         something_changed = true;
 451                         sign *= child_sign;
 452                 }
 453                 first++;
 454         }
 455
 456         // Now reorder the children
 457         first = symm.children.begin();
 458         switch (symm.type) {
 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));
 462                         break;
 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));
 466                         if (sign == 0)
 467                                 return 0;
 468                         break;
 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));
 472                         break;
 473                 default:
 474                         break;
 475         }
 476         return something_changed ? sign : std::numeric_limits<int>::max();
 477 }
 478
 479
 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)
 482 {
 483         // Need at least 2 objects for this operation
 484         unsigned num = last - first;
 485         if (num < 2)
 486                 return e;
 487
 488         // Transform object vector to a lst (for subs())
 489         lst orig_lst(first, last);
 490
 491         // Create index vectors for permutation
 492         unsigned *iv = new unsigned[num], *iv2;
 493         for (unsigned i=0; i<num; i++)
 494                 iv[i] = i;
 495         iv2 = (asymmetric ? new unsigned[num] : NULL);
 496
 497         // Loop over all permutations (the first permutation, which is the
 498         // identity, is unrolled)
 499         exvector sum_v;
 500         sum_v.push_back(e);
 501         while (std::next_permutation(iv, iv + num)) {
 502                 lst new_lst;
 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);
 506                 if (asymmetric) {
 507                         memcpy(iv2, iv, num * sizeof(unsigned));
 508                         term *= permutation_sign(iv2, iv2 + num);
 509                 }
 510                 sum_v.push_back(term);
 511         }
 512         ex sum = (new add(sum_v))->setflag(status_flags::dynallocated);
 513
 514         delete[] iv;
 515         delete[] iv2;
 516
 517         return sum / factorial(numeric(num));
 518 }
 519
 520 ex symmetrize(const ex & e, exvector::const_iterator first, exvector::const_iterator last)
 521 {
 522         return symm(e, first, last, false);
 523 }
 524
 525 ex antisymmetrize(const ex & e, exvector::const_iterator first, exvector::const_iterator last)
 526 {
 527         return symm(e, first, last, true);
 528 }
 529
 530 ex symmetrize_cyclic(const ex & e, exvector::const_iterator first, exvector::const_iterator last)
 531 {
 532         // Need at least 2 objects for this operation
 533         unsigned num = last - first;
 534         if (num < 2)
 535                 return e;
 536
 537         // Transform object vector to a lst (for subs())
 538         lst orig_lst(first, last);
 539         lst new_lst = orig_lst;
 540
 541         // Loop over all cyclic permutations (the first permutation, which is
 542         // the identity, is unrolled)
 543         ex sum = e;
 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);
 548         }
 549         return sum / num;
 550 }
 551
 552 /** Symmetrize expression over a list of objects (symbols, indices). */
 553 ex ex::symmetrize(const lst & l) const
 554 {
 555         exvector v(l.begin(), l.end());
 556         return symm(*this, v.begin(), v.end(), false);
 557 }
 558
 559 /** Antisymmetrize expression over a list of objects (symbols, indices). */
 560 ex ex::antisymmetrize(const lst & l) const
 561 {
 562         exvector v(l.begin(), l.end());
 563         return symm(*this, v.begin(), v.end(), true);
 564 }
 565
 566 /** Symmetrize expression by cyclic permutation over a list of objects
 567  *  (symbols, indices). */
 568 ex ex::symmetrize_cyclic(const lst & l) const
 569 {
 570         exvector v(l.begin(), l.end());
 571         return GiNaC::symmetrize_cyclic(*this, v.begin(), v.end());
 572 }
 573
 574 } // namespace GiNaC