ginac/symmetry.cpp

   1 /** @file symmetry.cpp
   2  *
   3  *  Implementation of GiNaC's symmetry definitions. */
   4
   5 /*
   6  *  GiNaC Copyright (C) 1999-2015 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_nonsymmetric() const
 274 {
 275         if (type == antisymmetric || 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_nonsymmetric())
 280                         return true;
 281
 282         return false;
 283 }
 284
 285 bool symmetry::has_cyclic() const
 286 {
 287         if (type == cyclic)
 288                 return true;
 289
 290         for (exvector::const_iterator i=children.begin(); i!=children.end(); ++i)
 291                 if (ex_to<symmetry>(*i).has_cyclic())
 292                         return true;
 293
 294         return false;
 295 }
 296
 297 symmetry &symmetry::add(const symmetry &c)
 298 {
 299         // All children must have the same number of indices
 300         if (type != none && !children.empty()) {
 301                 GINAC_ASSERT(is_exactly_a<symmetry>(children[0]));
 302                 if (ex_to<symmetry>(children[0]).indices.size() != c.indices.size())
 303                         throw (std::logic_error("symmetry:add(): children must have same number of indices"));
 304         }
 305
 306         // Compute union of indices and check whether the two sets are disjoint
 307         std::set<unsigned> un;
 308         set_union(indices.begin(), indices.end(), c.indices.begin(), c.indices.end(), inserter(un, un.begin()));
 309         if (un.size() != indices.size() + c.indices.size())
 310                 throw (std::logic_error("symmetry::add(): the same index appears in more than one child"));
 311
 312         // Set new index set
 313         indices.swap(un);
 314
 315         // Add child node
 316         children.push_back(c);
 317         return *this;
 318 }
 319
 320 void symmetry::validate(unsigned n)
 321 {
 322         if (indices.upper_bound(n - 1) != indices.end())
 323                 throw (std::range_error("symmetry::verify(): index values are out of range"));
 324         if (type != none && indices.empty()) {
 325                 for (unsigned i=0; i<n; i++)
 326                         add(i);
 327         }
 328 }
 329
 330 //////////
 331 // global functions
 332 //////////
 333
 334 static const symmetry & index0()
 335 {
 336         static ex s = (new symmetry(0))->setflag(status_flags::dynallocated);
 337         return ex_to<symmetry>(s);
 338 }
 339
 340 static const symmetry & index1()
 341 {
 342         static ex s = (new symmetry(1))->setflag(status_flags::dynallocated);
 343         return ex_to<symmetry>(s);
 344 }
 345
 346 static const symmetry & index2()
 347 {
 348         static ex s = (new symmetry(2))->setflag(status_flags::dynallocated);
 349         return ex_to<symmetry>(s);
 350 }
 351
 352 static const symmetry & index3()
 353 {
 354         static ex s = (new symmetry(3))->setflag(status_flags::dynallocated);
 355         return ex_to<symmetry>(s);
 356 }
 357
 358 const symmetry & not_symmetric()
 359 {
 360         static ex s = (new symmetry)->setflag(status_flags::dynallocated);
 361         return ex_to<symmetry>(s);
 362 }
 363
 364 const symmetry & symmetric2()
 365 {
 366         static ex s = (new symmetry(symmetry::symmetric, index0(), index1()))->setflag(status_flags::dynallocated);
 367         return ex_to<symmetry>(s);
 368 }
 369
 370 const symmetry & symmetric3()
 371 {
 372         static ex s = (new symmetry(symmetry::symmetric, index0(), index1()))->add(index2()).setflag(status_flags::dynallocated);
 373         return ex_to<symmetry>(s);
 374 }
 375
 376 const symmetry & symmetric4()
 377 {
 378         static ex s = (new symmetry(symmetry::symmetric, index0(), index1()))->add(index2()).add(index3()).setflag(status_flags::dynallocated);
 379         return ex_to<symmetry>(s);
 380 }
 381
 382 const symmetry & antisymmetric2()
 383 {
 384         static ex s = (new symmetry(symmetry::antisymmetric, index0(), index1()))->setflag(status_flags::dynallocated);
 385         return ex_to<symmetry>(s);
 386 }
 387
 388 const symmetry & antisymmetric3()
 389 {
 390         static ex s = (new symmetry(symmetry::antisymmetric, index0(), index1()))->add(index2()).setflag(status_flags::dynallocated);
 391         return ex_to<symmetry>(s);
 392 }
 393
 394 const symmetry & antisymmetric4()
 395 {
 396         static ex s = (new symmetry(symmetry::antisymmetric, index0(), index1()))->add(index2()).add(index3()).setflag(status_flags::dynallocated);
 397         return ex_to<symmetry>(s);
 398 }
 399
 400 class sy_is_less : public std::binary_function<ex, ex, bool> {
 401         exvector::iterator v;
 402
 403 public:
 404         sy_is_less(exvector::iterator v_) : v(v_) {}
 405
 406         bool operator() (const ex &lh, const ex &rh) const
 407         {
 408                 GINAC_ASSERT(is_exactly_a<symmetry>(lh));
 409                 GINAC_ASSERT(is_exactly_a<symmetry>(rh));
 410                 GINAC_ASSERT(ex_to<symmetry>(lh).indices.size() == ex_to<symmetry>(rh).indices.size());
 411                 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();
 412                 while (ait != aitend) {
 413                         int cmpval = v[*ait].compare(v[*bit]);
 414                         if (cmpval < 0)
 415                                 return true;
 416                         else if (cmpval > 0)
 417                                 return false;
 418                         ++ait; ++bit;
 419                 }
 420                 return false;
 421         }
 422 };
 423
 424 class sy_swap : public std::binary_function<ex, ex, void> {
 425         exvector::iterator v;
 426
 427 public:
 428         bool &swapped;
 429
 430         sy_swap(exvector::iterator v_, bool &s) : v(v_), swapped(s) {}
 431
 432         void operator() (const ex &lh, const ex &rh)
 433         {
 434                 GINAC_ASSERT(is_exactly_a<symmetry>(lh));
 435                 GINAC_ASSERT(is_exactly_a<symmetry>(rh));
 436                 GINAC_ASSERT(ex_to<symmetry>(lh).indices.size() == ex_to<symmetry>(rh).indices.size());
 437                 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();
 438                 while (ait != aitend) {
 439                         v[*ait].swap(v[*bit]);
 440                         ++ait; ++bit;
 441                 }
 442                 swapped = true;
 443         }
 444 };
 445
 446 int canonicalize(exvector::iterator v, const symmetry &symm)
 447 {
 448         // Less than two elements? Then do nothing
 449         if (symm.indices.size() < 2)
 450                 return std::numeric_limits<int>::max();
 451
 452         // Canonicalize children first
 453         bool something_changed = false;
 454         int sign = 1;
 455         exvector::const_iterator first = symm.children.begin(), last = symm.children.end();
 456         while (first != last) {
 457                 GINAC_ASSERT(is_exactly_a<symmetry>(*first));
 458                 int child_sign = canonicalize(v, ex_to<symmetry>(*first));
 459                 if (child_sign == 0)
 460                         return 0;
 461                 if (child_sign != std::numeric_limits<int>::max()) {
 462                         something_changed = true;
 463                         sign *= child_sign;
 464                 }
 465                 first++;
 466         }
 467
 468         // Now reorder the children
 469         first = symm.children.begin();
 470         switch (symm.type) {
 471                 case symmetry::symmetric:
 472                         // Sort the children in ascending order
 473                         shaker_sort(first, last, sy_is_less(v), sy_swap(v, something_changed));
 474                         break;
 475                 case symmetry::antisymmetric:
 476                         // Sort the children in ascending order, keeping track of the signum
 477                         sign *= permutation_sign(first, last, sy_is_less(v), sy_swap(v, something_changed));
 478                         if (sign == 0)
 479                                 return 0;
 480                         break;
 481                 case symmetry::cyclic:
 482                         // Permute the smallest child to the front
 483                         cyclic_permutation(first, last, min_element(first, last, sy_is_less(v)), sy_swap(v, something_changed));
 484                         break;
 485                 default:
 486                         break;
 487         }
 488         return something_changed ? sign : std::numeric_limits<int>::max();
 489 }
 490
 491
 492 // Symmetrize/antisymmetrize over a vector of objects
 493 static ex symm(const ex & e, exvector::const_iterator first, exvector::const_iterator last, bool asymmetric)
 494 {
 495         // Need at least 2 objects for this operation
 496         unsigned num = last - first;
 497         if (num < 2)
 498                 return e;
 499
 500         // Transform object vector to a lst (for subs())
 501         lst orig_lst(first, last);
 502
 503         // Create index vectors for permutation
 504         unsigned *iv = new unsigned[num], *iv2;
 505         for (unsigned i=0; i<num; i++)
 506                 iv[i] = i;
 507         iv2 = (asymmetric ? new unsigned[num] : nullptr);
 508
 509         // Loop over all permutations (the first permutation, which is the
 510         // identity, is unrolled)
 511         exvector sum_v;
 512         sum_v.push_back(e);
 513         while (std::next_permutation(iv, iv + num)) {
 514                 lst new_lst;
 515                 for (unsigned i=0; i<num; i++)
 516                         new_lst.append(orig_lst.op(iv[i]));
 517                 ex term = e.subs(orig_lst, new_lst, subs_options::no_pattern|subs_options::no_index_renaming);
 518                 if (asymmetric) {
 519                         memcpy(iv2, iv, num * sizeof(unsigned));
 520                         term *= permutation_sign(iv2, iv2 + num);
 521                 }
 522                 sum_v.push_back(term);
 523         }
 524         ex sum = (new add(sum_v))->setflag(status_flags::dynallocated);
 525
 526         delete[] iv;
 527         delete[] iv2;
 528
 529         return sum / factorial(numeric(num));
 530 }
 531
 532 ex symmetrize(const ex & e, exvector::const_iterator first, exvector::const_iterator last)
 533 {
 534         return symm(e, first, last, false);
 535 }
 536
 537 ex antisymmetrize(const ex & e, exvector::const_iterator first, exvector::const_iterator last)
 538 {
 539         return symm(e, first, last, true);
 540 }
 541
 542 ex symmetrize_cyclic(const ex & e, exvector::const_iterator first, exvector::const_iterator last)
 543 {
 544         // Need at least 2 objects for this operation
 545         unsigned num = last - first;
 546         if (num < 2)
 547                 return e;
 548
 549         // Transform object vector to a lst (for subs())
 550         lst orig_lst(first, last);
 551         lst new_lst = orig_lst;
 552
 553         // Loop over all cyclic permutations (the first permutation, which is
 554         // the identity, is unrolled)
 555         ex sum = e;
 556         for (unsigned i=0; i<num-1; i++) {
 557                 ex perm = new_lst.op(0);
 558                 new_lst.remove_first().append(perm);
 559                 sum += e.subs(orig_lst, new_lst, subs_options::no_pattern|subs_options::no_index_renaming);
 560         }
 561         return sum / num;
 562 }
 563
 564 /** Symmetrize expression over a list of objects (symbols, indices). */
 565 ex ex::symmetrize(const lst & l) const
 566 {
 567         exvector v(l.begin(), l.end());
 568         return symm(*this, v.begin(), v.end(), false);
 569 }
 570
 571 /** Antisymmetrize expression over a list of objects (symbols, indices). */
 572 ex ex::antisymmetrize(const lst & l) const
 573 {
 574         exvector v(l.begin(), l.end());
 575         return symm(*this, v.begin(), v.end(), true);
 576 }
 577
 578 /** Symmetrize expression by cyclic permutation over a list of objects
 579  *  (symbols, indices). */
 580 ex ex::symmetrize_cyclic(const lst & l) const
 581 {
 582         exvector v(l.begin(), l.end());
 583         return GiNaC::symmetrize_cyclic(*this, v.begin(), v.end());
 584 }
 585
 586 } // namespace GiNaC