GiNaC  1.8.0
indexed.cpp
Go to the documentation of this file.
1 
5 /*
6  * GiNaC Copyright (C) 1999-2020 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 "indexed.h"
24 #include "idx.h"
25 #include "add.h"
26 #include "mul.h"
27 #include "ncmul.h"
28 #include "power.h"
29 #include "relational.h"
30 #include "symmetry.h"
31 #include "operators.h"
32 #include "lst.h"
33 #include "archive.h"
34 #include "symbol.h"
35 #include "utils.h"
36 #include "integral.h"
37 #include "matrix.h"
38 #include "inifcns.h"
39 
40 #include <iostream>
41 #include <limits>
42 #include <sstream>
43 #include <stdexcept>
44 
45 namespace GiNaC {
46 
49  print_func<print_latex>(&indexed::do_print_latex).
50  print_func<print_tree>(&indexed::do_print_tree))
51 
52 // default constructor
55 
56 indexed::indexed() : symtree(not_symmetric())
57 {
58 }
59 
61 // other constructors
63 
64 indexed::indexed(const ex & b) : inherited{b}, symtree(not_symmetric())
65 {
66  validate();
67 }
68 
69 indexed::indexed(const ex & b, const ex & i1) : inherited{b, i1}, symtree(not_symmetric())
70 {
71  validate();
72 }
73 
74 indexed::indexed(const ex & b, const ex & i1, const ex & i2) : inherited{b, i1, i2}, symtree(not_symmetric())
75 {
76  validate();
77 }
78 
79 indexed::indexed(const ex & b, const ex & i1, const ex & i2, const ex & i3) : inherited{b, i1, i2, i3}, symtree(not_symmetric())
80 {
81  validate();
82 }
83 
84 indexed::indexed(const ex & b, const ex & i1, const ex & i2, const ex & i3, const ex & i4) : inherited{b, i1, i2, i3, i4}, symtree(not_symmetric())
85 {
86  validate();
87 }
88 
89 indexed::indexed(const ex & b, const symmetry & symm, const ex & i1, const ex & i2) : inherited{b, i1, i2}, symtree(symm)
90 {
91  validate();
92 }
93 
94 indexed::indexed(const ex & b, const symmetry & symm, const ex & i1, const ex & i2, const ex & i3) : inherited{b, i1, i2, i3}, symtree(symm)
95 {
96  validate();
97 }
98 
99 indexed::indexed(const ex & b, const symmetry & symm, const ex & i1, const ex & i2, const ex & i3, const ex & i4) : inherited{b, i1, i2, i3, i4}, symtree(symm)
100 {
101  validate();
102 }
103 
104 indexed::indexed(const ex & b, const exvector & v) : inherited{b}, symtree(not_symmetric())
105 {
106  seq.insert(seq.end(), v.begin(), v.end());
107  validate();
108 }
109 
110 indexed::indexed(const ex & b, const symmetry & symm, const exvector & v) : inherited{b}, symtree(symm)
111 {
112  seq.insert(seq.end(), v.begin(), v.end());
113  validate();
114 }
115 
116 indexed::indexed(const symmetry & symm, const exprseq & es) : inherited(es), symtree(symm)
117 {
118 }
119 
120 indexed::indexed(const symmetry & symm, const exvector & v) : inherited(v), symtree(symm)
121 {
122 }
123 
124 indexed::indexed(const symmetry & symm, exvector && v) : inherited(std::move(v)), symtree(symm)
125 {
126 }
127 
129 // archiving
131 
132 void indexed::read_archive(const archive_node &n, lst &sym_lst)
133 {
134  inherited::read_archive(n, sym_lst);
135  if (!n.find_ex("symmetry", symtree, sym_lst)) {
136  // GiNaC versions <= 0.9.0 had an unsigned "symmetry" property
137  unsigned symm = 0;
138  n.find_unsigned("symmetry", symm);
139  switch (symm) {
140  case 1:
141  symtree = sy_symm();
142  break;
143  case 2:
144  symtree = sy_anti();
145  break;
146  default:
148  break;
149  }
150  const_cast<symmetry &>(ex_to<symmetry>(symtree)).validate(seq.size() - 1);
151  }
152 }
154 
156 {
157  inherited::archive(n);
158  n.add_ex("symmetry", symtree);
159 }
160 
162 // functions overriding virtual functions from base classes
164 
165 void indexed::printindices(const print_context & c, unsigned level) const
166 {
167  if (seq.size() > 1) {
168 
169  auto it = seq.begin() + 1, itend = seq.end();
170 
171  if (is_a<print_latex>(c)) {
172 
173  // TeX output: group by variance
174  bool first = true;
175  bool covariant = true;
176 
177  while (it != itend) {
178  bool cur_covariant = (is_a<varidx>(*it) ? ex_to<varidx>(*it).is_covariant() : true);
179  if (first || cur_covariant != covariant) { // Variance changed
180  // The empty {} prevents indices from ending up on top of each other
181  if (!first)
182  c.s << "}{}";
183  covariant = cur_covariant;
184  if (covariant)
185  c.s << "_{";
186  else
187  c.s << "^{";
188  }
189  it->print(c, level);
190  c.s << " ";
191  first = false;
192  it++;
193  }
194  c.s << "}";
195 
196  } else {
197 
198  // Ordinary output
199  while (it != itend) {
200  it->print(c, level);
201  it++;
202  }
203  }
204  }
205 }
206 
207 void indexed::print_indexed(const print_context & c, const char *openbrace, const char *closebrace, unsigned level) const
208 {
209  if (precedence() <= level)
210  c.s << openbrace << '(';
211  c.s << openbrace;
212  seq[0].print(c, precedence());
213  c.s << closebrace;
214  printindices(c, level);
215  if (precedence() <= level)
216  c.s << ')' << closebrace;
217 }
218 
219 void indexed::do_print(const print_context & c, unsigned level) const
220 {
221  print_indexed(c, "", "", level);
222 }
223 
224 void indexed::do_print_latex(const print_latex & c, unsigned level) const
225 {
226  print_indexed(c, "{", "}", level);
227 }
228 
229 void indexed::do_print_tree(const print_tree & c, unsigned level) const
230 {
231  c.s << std::string(level, ' ') << class_name() << " @" << this
232  << std::hex << ", hash=0x" << hashvalue << ", flags=0x" << flags << std::dec
233  << ", " << seq.size()-1 << " indices"
234  << ", symmetry=" << symtree << std::endl;
235  seq[0].print(c, level + c.delta_indent);
236  printindices(c, level + c.delta_indent);
237 }
238 
239 bool indexed::info(unsigned inf) const
240 {
241  if (inf == info_flags::indexed) return true;
242  if (inf == info_flags::has_indices) return seq.size() > 1;
243  return inherited::info(inf);
244 }
245 
246 bool indexed::all_index_values_are(unsigned inf) const
247 {
248  // No indices? Then no property can be fulfilled
249  if (seq.size() < 2)
250  return false;
251 
252  // Check all indices
253  return find_if(seq.begin() + 1, seq.end(),
254  [inf](const ex & e) { return !(ex_to<idx>(e).get_value().info(inf)); }) == seq.end();
255 }
256 
257 int indexed::compare_same_type(const basic & other) const
258 {
259  GINAC_ASSERT(is_a<indexed>(other));
260  return inherited::compare_same_type(other);
261 }
262 
264 {
265  const ex &base = seq[0];
266 
267  // If the base object is 0, the whole object is 0
268  if (base.is_zero())
269  return _ex0;
270 
271  // If the base object is a product, pull out the numeric factor
272  if (is_exactly_a<mul>(base) && is_exactly_a<numeric>(base.op(base.nops() - 1))) {
273  exvector v(seq);
274  ex f = ex_to<numeric>(base.op(base.nops() - 1));
275  v[0] = seq[0] / f;
276  return f * thiscontainer(v);
277  }
278 
279  if((typeid(*this) == typeid(indexed)) && seq.size()==1)
280  return base;
281 
282  // Canonicalize indices according to the symmetry properties
283  if (seq.size() > 2) {
284  exvector v = seq;
285  GINAC_ASSERT(is_exactly_a<symmetry>(symtree));
286  int sig = canonicalize(v.begin() + 1, ex_to<symmetry>(symtree));
287  if (sig != std::numeric_limits<int>::max()) {
288  // Something has changed while sorting indices, more evaluations later
289  if (sig == 0)
290  return _ex0;
291  return ex(sig) * thiscontainer(v);
292  }
293  }
294 
295  // Let the class of the base object perform additional evaluations
296  return ex_to<basic>(base).eval_indexed(*this);
297 }
298 
300 {
301  if(op(0).info(info_flags::real))
302  return *this;
303  return real_part_function(*this).hold();
304 }
305 
307 {
308  if(op(0).info(info_flags::real))
309  return 0;
310  return imag_part_function(*this).hold();
311 }
312 
314 {
315  return indexed(ex_to<symmetry>(symtree), v);
316 }
317 
319 {
320  return indexed(ex_to<symmetry>(symtree), std::move(v));
321 }
322 
323 unsigned indexed::return_type() const
324 {
325  if(is_a<matrix>(op(0)))
327  else
328  return op(0).return_type();
329 }
330 
331 ex indexed::expand(unsigned options) const
332 {
333  GINAC_ASSERT(seq.size() > 0);
334 
336  ex newbase = seq[0].expand(options);
337  if (is_exactly_a<add>(newbase)) {
338  ex sum = _ex0;
339  for (size_t i=0; i<newbase.nops(); i++) {
340  exvector s = seq;
341  s[0] = newbase.op(i);
342  sum += thiscontainer(s).expand(options);
343  }
344  return sum;
345  }
346  if (!are_ex_trivially_equal(newbase, seq[0])) {
347  exvector s = seq;
348  s[0] = newbase;
349  return ex_to<indexed>(thiscontainer(s)).inherited::expand(options);
350  }
351  }
352  return inherited::expand(options);
353 }
354 
356 // virtual functions which can be overridden by derived classes
358 
359 // none
360 
362 // non-virtual functions in this class
364 
368 void indexed::validate() const
369 {
370  GINAC_ASSERT(seq.size() > 0);
371  auto it = seq.begin() + 1, itend = seq.end();
372  while (it != itend) {
373  if (!is_a<idx>(*it))
374  throw(std::invalid_argument("indices of indexed object must be of type idx"));
375  it++;
376  }
377 
378  if (!symtree.is_zero()) {
379  if (!is_exactly_a<symmetry>(symtree))
380  throw(std::invalid_argument("symmetry of indexed object must be of type symmetry"));
381  const_cast<symmetry &>(ex_to<symmetry>(symtree)).validate(seq.size() - 1);
382  }
383 }
384 
388 ex indexed::derivative(const symbol & s) const
389 {
390  return _ex0;
391 }
392 
394 // global functions
396 
398  bool operator() (const ex &lh, const ex &rh) const
399  {
400  if (lh.is_equal(rh))
401  return true;
402  else
403  try {
404  // Replacing the dimension might cause an error (e.g. with
405  // index classes that only work in a fixed number of dimensions)
406  return lh.is_equal(ex_to<idx>(rh).replace_dim(ex_to<idx>(lh).get_dim()));
407  } catch (...) {
408  return false;
409  }
410  }
411 };
412 
414 static bool indices_consistent(const exvector & v1, const exvector & v2)
415 {
416  // Number of indices must be the same
417  if (v1.size() != v2.size())
418  return false;
419 
420  return equal(v1.begin(), v1.end(), v2.begin(), idx_is_equal_ignore_dim());
421 }
422 
424 {
425  GINAC_ASSERT(seq.size() >= 1);
426  return exvector(seq.begin() + 1, seq.end());
427 }
428 
430 {
431  exvector free_indices, dummy_indices;
432  find_free_and_dummy(seq.begin() + 1, seq.end(), free_indices, dummy_indices);
433  return dummy_indices;
434 }
435 
437 {
438  exvector indices = get_free_indices();
439  exvector other_indices = other.get_free_indices();
440  indices.insert(indices.end(), other_indices.begin(), other_indices.end());
441  exvector dummy_indices;
442  find_dummy_indices(indices, dummy_indices);
443  return dummy_indices;
444 }
445 
446 bool indexed::has_dummy_index_for(const ex & i) const
447 {
448  auto it = seq.begin() + 1, itend = seq.end();
449  while (it != itend) {
450  if (is_dummy_pair(*it, i))
451  return true;
452  it++;
453  }
454  return false;
455 }
456 
458 {
459  exvector free_indices, dummy_indices;
460  find_free_and_dummy(seq.begin() + 1, seq.end(), free_indices, dummy_indices);
461  return free_indices;
462 }
463 
465 {
466  exvector free_indices;
467  for (size_t i=0; i<nops(); i++) {
468  if (i == 0)
469  free_indices = op(i).get_free_indices();
470  else {
471  exvector free_indices_of_term = op(i).get_free_indices();
472  if (!indices_consistent(free_indices, free_indices_of_term))
473  throw (std::runtime_error("add::get_free_indices: inconsistent indices in sum"));
474  }
475  }
476  return free_indices;
477 }
478 
480 {
481  // Concatenate free indices of all factors
482  exvector un;
483  for (size_t i=0; i<nops(); i++) {
484  exvector free_indices_of_factor = op(i).get_free_indices();
485  un.insert(un.end(), free_indices_of_factor.begin(), free_indices_of_factor.end());
486  }
487 
488  // And remove the dummy indices
489  exvector free_indices, dummy_indices;
490  find_free_and_dummy(un, free_indices, dummy_indices);
491  return free_indices;
492 }
493 
495 {
496  // Concatenate free indices of all factors
497  exvector un;
498  for (size_t i=0; i<nops(); i++) {
499  exvector free_indices_of_factor = op(i).get_free_indices();
500  un.insert(un.end(), free_indices_of_factor.begin(), free_indices_of_factor.end());
501  }
502 
503  // And remove the dummy indices
504  exvector free_indices, dummy_indices;
505  find_free_and_dummy(un, free_indices, dummy_indices);
506  return free_indices;
507 }
508 
510  bool operator()(const ex & e)
511  {
512  return is_dummy_pair(e, e);
513  }
514 };
515 
517 {
518  if (a.get_free_indices().size() || b.get_free_indices().size())
519  throw (std::runtime_error("integral::get_free_indices: boundary values should not have free indices"));
520  return f.get_free_indices();
521 }
522 
523 template<class T> size_t number_of_type(const exvector&v)
524 {
525  size_t number = 0;
526  for (auto & it : v)
527  if (is_exactly_a<T>(it))
528  ++number;
529  return number;
530 }
531 
540 template<class T> static ex rename_dummy_indices(const ex & e, exvector & global_dummy_indices, exvector & local_dummy_indices)
541 {
542  size_t global_size = number_of_type<T>(global_dummy_indices),
543  local_size = number_of_type<T>(local_dummy_indices);
544 
545  // Any local dummy indices at all?
546  if (local_size == 0)
547  return e;
548 
549  if (global_size < local_size) {
550 
551  // More local indices than we encountered before, add the new ones
552  // to the global set
553  size_t old_global_size = global_size;
554  int remaining = local_size - global_size;
555  auto it = local_dummy_indices.begin(), itend = local_dummy_indices.end();
556  while (it != itend && remaining > 0) {
557  if (is_exactly_a<T>(*it) &&
558  find_if(global_dummy_indices.begin(), global_dummy_indices.end(),
559  [it](const ex &lh) { return idx_is_equal_ignore_dim()(lh, *it); }) == global_dummy_indices.end()) {
560  global_dummy_indices.push_back(*it);
561  global_size++;
562  remaining--;
563  }
564  it++;
565  }
566 
567  // If this is the first set of local indices, do nothing
568  if (old_global_size == 0)
569  return e;
570  }
571  GINAC_ASSERT(local_size <= global_size);
572 
573  // Construct vectors of index symbols
574  exvector local_syms, global_syms;
575  local_syms.reserve(local_size);
576  global_syms.reserve(local_size);
577  for (size_t i=0; local_syms.size()!=local_size; i++)
578  if(is_exactly_a<T>(local_dummy_indices[i]))
579  local_syms.push_back(local_dummy_indices[i].op(0));
580  shaker_sort(local_syms.begin(), local_syms.end(), ex_is_less(), ex_swap());
581  for (size_t i=0; global_syms.size()!=local_size; i++) // don't use more global symbols than necessary
582  if(is_exactly_a<T>(global_dummy_indices[i]))
583  global_syms.push_back(global_dummy_indices[i].op(0));
584  shaker_sort(global_syms.begin(), global_syms.end(), ex_is_less(), ex_swap());
585 
586  // Remove common indices
587  exvector local_uniq, global_uniq;
588  set_difference(local_syms.begin(), local_syms.end(), global_syms.begin(), global_syms.end(), std::back_insert_iterator<exvector>(local_uniq), ex_is_less());
589  set_difference(global_syms.begin(), global_syms.end(), local_syms.begin(), local_syms.end(), std::back_insert_iterator<exvector>(global_uniq), ex_is_less());
590 
591  // Replace remaining non-common local index symbols by global ones
592  if (local_uniq.empty())
593  return e;
594  else {
595  while (global_uniq.size() > local_uniq.size())
596  global_uniq.pop_back();
597  return e.subs(lst(local_uniq.begin(), local_uniq.end()), lst(global_uniq.begin(), global_uniq.end()), subs_options::no_pattern);
598  }
599 }
600 
602 static void find_variant_indices(const exvector & v, exvector & variant_indices)
603 {
604  exvector::const_iterator it1, itend;
605  for (it1 = v.begin(), itend = v.end(); it1 != itend; ++it1) {
606  if (is_exactly_a<varidx>(*it1))
607  variant_indices.push_back(*it1);
608  }
609 }
610 
618 bool reposition_dummy_indices(ex & e, exvector & variant_dummy_indices, exvector & moved_indices)
619 {
620  bool something_changed = false;
621 
622  // Find dummy symbols that occur twice in the same indexed object.
623  exvector local_var_dummies;
624  local_var_dummies.reserve(e.nops()/2);
625  for (size_t i=1; i<e.nops(); ++i) {
626  if (!is_a<varidx>(e.op(i)))
627  continue;
628  for (size_t j=i+1; j<e.nops(); ++j) {
629  if (is_dummy_pair(e.op(i), e.op(j))) {
630  local_var_dummies.push_back(e.op(i));
631  for (auto k = variant_dummy_indices.begin(); k!=variant_dummy_indices.end(); ++k) {
632  if (e.op(i).op(0) == k->op(0)) {
633  variant_dummy_indices.erase(k);
634  break;
635  }
636  }
637  break;
638  }
639  }
640  }
641 
642  // In the case where a dummy symbol occurs twice in the same indexed object
643  // we try all possibilities of raising/lowering and keep the least one in
644  // the sense of ex_is_less.
645  ex optimal_e = e;
646  size_t numpossibs = 1 << local_var_dummies.size();
647  for (size_t i=0; i<numpossibs; ++i) {
648  ex try_e = e;
649  for (size_t j=0; j<local_var_dummies.size(); ++j) {
650  exmap m;
651  if (1<<j & i) {
652  ex curr_idx = local_var_dummies[j];
653  ex curr_toggle = ex_to<varidx>(curr_idx).toggle_variance();
654  m[curr_idx] = curr_toggle;
655  m[curr_toggle] = curr_idx;
656  }
657  try_e = e.subs(m, subs_options::no_pattern);
658  }
659  if(ex_is_less()(try_e, optimal_e))
660  { optimal_e = try_e;
661  something_changed = true;
662  }
663  }
664  e = optimal_e;
665 
666  if (!is_a<indexed>(e))
667  return true;
668 
669  exvector seq = ex_to<indexed>(e).seq;
670 
671  // If a dummy index is encountered for the first time in the
672  // product, pull it up, otherwise, pull it down
673  for (auto it2 = seq.begin()+1, it2end = seq.end(); it2 != it2end; ++it2) {
674  if (!is_exactly_a<varidx>(*it2))
675  continue;
676 
677  exvector::iterator vit, vitend;
678  for (vit = variant_dummy_indices.begin(), vitend = variant_dummy_indices.end(); vit != vitend; ++vit) {
679  if (it2->op(0).is_equal(vit->op(0))) {
680  if (ex_to<varidx>(*it2).is_covariant()) {
681  /*
682  * N.B. we don't want to use
683  *
684  * e = e.subs(lst{
685  * *it2 == ex_to<varidx>(*it2).toggle_variance(),
686  * ex_to<varidx>(*it2).toggle_variance() == *it2
687  * }, subs_options::no_pattern);
688  *
689  * since this can trigger non-trivial repositioning of indices,
690  * e.g. due to non-trivial symmetry properties of e, thus
691  * invalidating iterators
692  */
693  *it2 = ex_to<varidx>(*it2).toggle_variance();
694  something_changed = true;
695  }
696  moved_indices.push_back(*vit);
697  variant_dummy_indices.erase(vit);
698  goto next_index;
699  }
700  }
701 
702  for (vit = moved_indices.begin(), vitend = moved_indices.end(); vit != vitend; ++vit) {
703  if (it2->op(0).is_equal(vit->op(0))) {
704  if (ex_to<varidx>(*it2).is_contravariant()) {
705  *it2 = ex_to<varidx>(*it2).toggle_variance();
706  something_changed = true;
707  }
708  goto next_index;
709  }
710  }
711 
712 next_index: ;
713  }
714 
715  if (something_changed)
716  e = ex_to<indexed>(e).thiscontainer(seq);
717 
718  return something_changed;
719 }
720 
721 /* Ordering that only compares the base expressions of indexed objects. */
723  bool operator() (const ex &lh, const ex &rh) const
724  {
725  return (is_a<indexed>(lh) ? lh.op(0) : lh).compare(is_a<indexed>(rh) ? rh.op(0) : rh) < 0;
726  }
727 };
728 
729 /* An auxiliary function used by simplify_indexed() and expand_dummy_sum()
730  * It returns an exvector of factors from the supplied product */
731 static void product_to_exvector(const ex & e, exvector & v, bool & non_commutative)
732 {
733  // Remember whether the product was commutative or noncommutative
734  // (because we chop it into factors and need to reassemble later)
735  non_commutative = is_exactly_a<ncmul>(e);
736 
737  // Collect factors in an exvector, store squares twice
738  v.reserve(e.nops() * 2);
739 
740  if (is_exactly_a<power>(e)) {
741  // We only get called for simple squares, split a^2 -> a*a
742  GINAC_ASSERT(e.op(1).is_equal(_ex2));
743  v.push_back(e.op(0));
744  v.push_back(e.op(0));
745  } else {
746  for (size_t i=0; i<e.nops(); i++) {
747  ex f = e.op(i);
748  if (is_exactly_a<power>(f) && f.op(1).is_equal(_ex2)) {
749  v.push_back(f.op(0));
750  v.push_back(f.op(0));
751  } else if (is_exactly_a<ncmul>(f)) {
752  // Noncommutative factor found, split it as well
753  non_commutative = true; // everything becomes noncommutative, ncmul will sort out the commutative factors later
754  for (size_t j=0; j<f.nops(); j++)
755  v.push_back(f.op(j));
756  } else
757  v.push_back(f);
758  }
759  }
760 }
761 
762 template<class T> ex idx_symmetrization(const ex& r,const exvector& local_dummy_indices)
763 { exvector dummy_syms;
764  dummy_syms.reserve(r.nops());
765  for (auto & it : local_dummy_indices)
766  if(is_exactly_a<T>(it))
767  dummy_syms.push_back(it.op(0));
768  if(dummy_syms.size() < 2)
769  return r;
770  ex q=symmetrize(r, dummy_syms);
771  return q;
772 }
773 
774 // Forward declaration needed in absence of friend injection, C.f. [namespace.memdef]:
775 ex simplify_indexed(const ex & e, exvector & free_indices, exvector & dummy_indices, const scalar_products & sp);
776 
779 ex simplify_indexed_product(const ex & e, exvector & free_indices, exvector & dummy_indices, const scalar_products & sp)
780 {
781  // Collect factors in an exvector
782  exvector v;
783 
784  // Remember whether the product was commutative or noncommutative
785  // (because we chop it into factors and need to reassemble later)
786  bool non_commutative;
787  product_to_exvector(e, v, non_commutative);
788 
789  // Perform contractions
790  bool something_changed = false;
791  bool has_nonsymmetric = false;
792  GINAC_ASSERT(v.size() > 1);
793  exvector::iterator it1, itend = v.end(), next_to_last = itend - 1;
794  for (it1 = v.begin(); it1 != next_to_last; it1++) {
795 
796 try_again:
797  if (!is_a<indexed>(*it1))
798  continue;
799 
800  bool first_noncommutative = (it1->return_type() != return_types::commutative);
801  bool first_nonsymmetric = ex_to<symmetry>(ex_to<indexed>(*it1).get_symmetry()).has_nonsymmetric();
802 
803  // Indexed factor found, get free indices and look for contraction
804  // candidates
805  exvector free1, dummy1;
806  find_free_and_dummy(ex_to<indexed>(*it1).seq.begin() + 1, ex_to<indexed>(*it1).seq.end(), free1, dummy1);
807 
808  exvector::iterator it2;
809  for (it2 = it1 + 1; it2 != itend; it2++) {
810 
811  if (!is_a<indexed>(*it2))
812  continue;
813 
814  bool second_noncommutative = (it2->return_type() != return_types::commutative);
815 
816  // Find free indices of second factor and merge them with free
817  // indices of first factor
818  exvector un;
819  find_free_and_dummy(ex_to<indexed>(*it2).seq.begin() + 1, ex_to<indexed>(*it2).seq.end(), un, dummy1);
820  un.insert(un.end(), free1.begin(), free1.end());
821 
822  // Check whether the two factors share dummy indices
823  exvector free, dummy;
824  find_free_and_dummy(un, free, dummy);
825  size_t num_dummies = dummy.size();
826  if (num_dummies == 0)
827  continue;
828 
829  // At least one dummy index, is it a defined scalar product?
830  bool contracted = false;
831  if (free.empty() && it1->nops()==2 && it2->nops()==2) {
832 
833  ex dim = minimal_dim(
834  ex_to<idx>(it1->op(1)).get_dim(),
835  ex_to<idx>(it2->op(1)).get_dim()
836  );
837 
838  // User-defined scalar product?
839  if (sp.is_defined(*it1, *it2, dim)) {
840 
841  // Yes, substitute it
842  *it1 = sp.evaluate(*it1, *it2, dim);
843  *it2 = _ex1;
844  goto contraction_done;
845  }
846  }
847 
848  // Try to contract the first one with the second one
849  contracted = ex_to<basic>(it1->op(0)).contract_with(it1, it2, v);
850  if (!contracted) {
851 
852  // That didn't work; maybe the second object knows how to
853  // contract itself with the first one
854  contracted = ex_to<basic>(it2->op(0)).contract_with(it2, it1, v);
855  }
856  if (contracted) {
857 contraction_done:
858  if (first_noncommutative || second_noncommutative
859  || is_exactly_a<add>(*it1) || is_exactly_a<add>(*it2)
860  || is_exactly_a<mul>(*it1) || is_exactly_a<mul>(*it2)
861  || is_exactly_a<ncmul>(*it1) || is_exactly_a<ncmul>(*it2)) {
862 
863  // One of the factors became a sum or product:
864  // re-expand expression and run again
865  // Non-commutative products are always re-expanded to give
866  // eval_ncmul() the chance to re-order and canonicalize
867  // the product
868  bool is_a_product = (is_exactly_a<mul>(*it1) || is_exactly_a<ncmul>(*it1)) &&
869  (is_exactly_a<mul>(*it2) || is_exactly_a<ncmul>(*it2));
870  ex r = (non_commutative ? ex(ncmul(std::move(v))) : ex(mul(std::move(v))));
871 
872  // If new expression is a product we can call this function again,
873  // otherwise we need to pass argument to simplify_indexed() to be expanded
874  if (is_a_product)
875  return simplify_indexed_product(r, free_indices, dummy_indices, sp);
876  else
877  return simplify_indexed(r, free_indices, dummy_indices, sp);
878  }
879 
880  // Both objects may have new indices now or they might
881  // even not be indexed objects any more, so we have to
882  // start over
883  something_changed = true;
884  goto try_again;
885  }
886  else if (!has_nonsymmetric &&
887  (first_nonsymmetric ||
888  ex_to<symmetry>(ex_to<indexed>(*it2).get_symmetry()).has_nonsymmetric())) {
889  has_nonsymmetric = true;
890  }
891  }
892  }
893 
894  // Find free indices (concatenate them all and call find_free_and_dummy())
895  // and all dummy indices that appear
896  exvector un, individual_dummy_indices;
897  for (it1 = v.begin(), itend = v.end(); it1 != itend; ++it1) {
898  exvector free_indices_of_factor;
899  if (is_a<indexed>(*it1)) {
900  exvector dummy_indices_of_factor;
901  find_free_and_dummy(ex_to<indexed>(*it1).seq.begin() + 1, ex_to<indexed>(*it1).seq.end(), free_indices_of_factor, dummy_indices_of_factor);
902  individual_dummy_indices.insert(individual_dummy_indices.end(), dummy_indices_of_factor.begin(), dummy_indices_of_factor.end());
903  } else
904  free_indices_of_factor = it1->get_free_indices();
905  un.insert(un.end(), free_indices_of_factor.begin(), free_indices_of_factor.end());
906  }
907  exvector local_dummy_indices;
908  find_free_and_dummy(un, free_indices, local_dummy_indices);
909  local_dummy_indices.insert(local_dummy_indices.end(), individual_dummy_indices.begin(), individual_dummy_indices.end());
910 
911  // Filter out the dummy indices with variance
912  exvector variant_dummy_indices;
913  find_variant_indices(local_dummy_indices, variant_dummy_indices);
914 
915  // Any indices with variance present at all?
916  if (!variant_dummy_indices.empty()) {
917 
918  // Yes, bring the product into a canonical order that only depends on
919  // the base expressions of indexed objects
920  if (!non_commutative)
921  std::sort(v.begin(), v.end(), ex_base_is_less());
922 
923  exvector moved_indices;
924 
925  // Iterate over all indexed objects in the product
926  for (it1 = v.begin(), itend = v.end(); it1 != itend; ++it1) {
927  if (!is_a<indexed>(*it1))
928  continue;
929 
930  if (reposition_dummy_indices(*it1, variant_dummy_indices, moved_indices))
931  something_changed = true;
932  }
933  }
934 
935  ex r;
936  if (something_changed)
937  r = non_commutative ? ex(ncmul(std::move(v))) : ex(mul(std::move(v)));
938  else
939  r = e;
940 
941  // The result should be symmetric with respect to exchange of dummy
942  // indices, so if the symmetrization vanishes, the whole expression is
943  // zero. This detects things like eps.i.j.k * p.j * p.k = 0.
944  if (has_nonsymmetric) {
945  ex q = idx_symmetrization<idx>(r, local_dummy_indices);
946  if (q.is_zero()) {
947  free_indices.clear();
948  return _ex0;
949  }
950  q = idx_symmetrization<varidx>(q, local_dummy_indices);
951  if (q.is_zero()) {
952  free_indices.clear();
953  return _ex0;
954  }
955  q = idx_symmetrization<spinidx>(q, local_dummy_indices);
956  if (q.is_zero()) {
957  free_indices.clear();
958  return _ex0;
959  }
960  }
961 
962  // Dummy index renaming
963  r = rename_dummy_indices<idx>(r, dummy_indices, local_dummy_indices);
964  r = rename_dummy_indices<varidx>(r, dummy_indices, local_dummy_indices);
965  r = rename_dummy_indices<spinidx>(r, dummy_indices, local_dummy_indices);
966 
967  // Product of indexed object with a scalar?
968  if (is_exactly_a<mul>(r) && r.nops() == 2
969  && is_exactly_a<numeric>(r.op(1)) && is_a<indexed>(r.op(0)))
970  return ex_to<basic>(r.op(0).op(0)).scalar_mul_indexed(r.op(0), ex_to<numeric>(r.op(1)));
971  else
972  return r;
973 }
974 
977 class terminfo {
978 public:
979  terminfo(const ex & orig_, const ex & symm_) : orig(orig_), symm(symm_) {}
980 
983 };
984 
986 public:
987  bool operator() (const terminfo & ti1, const terminfo & ti2) const
988  {
989  return (ti1.symm.compare(ti2.symm) < 0);
990  }
991 };
992 
995 class symminfo {
996 public:
997  symminfo() : num(0) {}
998 
999  symminfo(const ex & symmterm_, const ex & orig_, size_t num_) : orig(orig_), num(num_)
1000  {
1001  if (is_exactly_a<mul>(symmterm_) && is_exactly_a<numeric>(symmterm_.op(symmterm_.nops()-1))) {
1002  coeff = symmterm_.op(symmterm_.nops()-1);
1003  symmterm = symmterm_ / coeff;
1004  } else {
1005  coeff = 1;
1006  symmterm = symmterm_;
1007  }
1008  }
1009 
1013  size_t num;
1014 };
1015 
1017 public:
1018  bool operator() (const symminfo & si1, const symminfo & si2) const
1019  {
1020  return (si1.symmterm.compare(si2.symmterm) < 0);
1021  }
1022 };
1023 
1025 public:
1026  bool operator() (const symminfo & si1, const symminfo & si2) const
1027  {
1028  return (si1.orig.compare(si2.orig) < 0);
1029  }
1030 };
1031 
1032 bool hasindex(const ex &x, const ex &sym)
1033 {
1034  if(is_a<idx>(x) && x.op(0)==sym)
1035  return true;
1036  else
1037  for(size_t i=0; i<x.nops(); ++i)
1038  if(hasindex(x.op(i), sym))
1039  return true;
1040  return false;
1041 }
1042 
1044 ex simplify_indexed(const ex & e, exvector & free_indices, exvector & dummy_indices, const scalar_products & sp)
1045 {
1046  // Expand the expression
1047  ex e_expanded = e.expand();
1048 
1049  // Simplification of single indexed object: just find the free indices
1050  // and perform dummy index renaming/repositioning
1051  if (is_a<indexed>(e_expanded)) {
1052 
1053  // Find the dummy indices
1054  const indexed &i = ex_to<indexed>(e_expanded);
1055  exvector local_dummy_indices;
1056  find_free_and_dummy(i.seq.begin() + 1, i.seq.end(), free_indices, local_dummy_indices);
1057 
1058  // Filter out the dummy indices with variance
1059  exvector variant_dummy_indices;
1060  find_variant_indices(local_dummy_indices, variant_dummy_indices);
1061 
1062  // Any indices with variance present at all?
1063  if (!variant_dummy_indices.empty()) {
1064 
1065  // Yes, reposition them
1066  exvector moved_indices;
1067  reposition_dummy_indices(e_expanded, variant_dummy_indices, moved_indices);
1068  }
1069 
1070  // Rename the dummy indices
1071  e_expanded = rename_dummy_indices<idx>(e_expanded, dummy_indices, local_dummy_indices);
1072  e_expanded = rename_dummy_indices<varidx>(e_expanded, dummy_indices, local_dummy_indices);
1073  e_expanded = rename_dummy_indices<spinidx>(e_expanded, dummy_indices, local_dummy_indices);
1074  return e_expanded;
1075  }
1076 
1077  // Simplification of sum = sum of simplifications, check consistency of
1078  // free indices in each term
1079  if (is_exactly_a<add>(e_expanded)) {
1080  bool first = true;
1081  ex sum;
1082  free_indices.clear();
1083 
1084  for (size_t i=0; i<e_expanded.nops(); i++) {
1085  exvector free_indices_of_term;
1086  ex term = simplify_indexed(e_expanded.op(i), free_indices_of_term, dummy_indices, sp);
1087  if (!term.is_zero()) {
1088  if (first) {
1089  free_indices = free_indices_of_term;
1090  sum = term;
1091  first = false;
1092  } else {
1093  if (!indices_consistent(free_indices, free_indices_of_term)) {
1094  std::ostringstream s;
1095  s << "simplify_indexed: inconsistent indices in sum: ";
1096  s << exprseq(free_indices) << " vs. " << exprseq(free_indices_of_term);
1097  throw (std::runtime_error(s.str()));
1098  }
1099  if (is_a<indexed>(sum) && is_a<indexed>(term))
1100  sum = ex_to<basic>(sum.op(0)).add_indexed(sum, term);
1101  else
1102  sum += term;
1103  }
1104  }
1105  }
1106 
1107  // If the sum turns out to be zero, we are finished
1108  if (sum.is_zero()) {
1109  free_indices.clear();
1110  return sum;
1111  }
1112 
1113  // More than one term and more than one dummy index?
1114  size_t num_terms_orig = (is_exactly_a<add>(sum) ? sum.nops() : 1);
1115  if (num_terms_orig < 2 || dummy_indices.size() < 2)
1116  return sum;
1117 
1118  // Chop the sum into terms and symmetrize each one over the dummy
1119  // indices
1120  std::vector<terminfo> terms;
1121  for (size_t i=0; i<sum.nops(); i++) {
1122  const ex & term = sum.op(i);
1123  exvector dummy_indices_of_term;
1124  dummy_indices_of_term.reserve(dummy_indices.size());
1125  for (auto & i : dummy_indices)
1126  if (hasindex(term,i.op(0)))
1127  dummy_indices_of_term.push_back(i);
1128  ex term_symm = idx_symmetrization<idx>(term, dummy_indices_of_term);
1129  term_symm = idx_symmetrization<varidx>(term_symm, dummy_indices_of_term);
1130  term_symm = idx_symmetrization<spinidx>(term_symm, dummy_indices_of_term);
1131  if (term_symm.is_zero())
1132  continue;
1133  terms.push_back(terminfo(term, term_symm));
1134  }
1135 
1136  // Sort by symmetrized terms
1137  std::sort(terms.begin(), terms.end(), terminfo_is_less());
1138 
1139  // Combine equal symmetrized terms
1140  std::vector<terminfo> terms_pass2;
1141  for (std::vector<terminfo>::const_iterator i=terms.begin(); i!=terms.end(); ) {
1142  size_t num = 1;
1143  auto j = i + 1;
1144  while (j != terms.end() && j->symm == i->symm) {
1145  num++;
1146  j++;
1147  }
1148  terms_pass2.push_back(terminfo(i->orig * num, i->symm * num));
1149  i = j;
1150  }
1151 
1152  // If there is only one term left, we are finished
1153  if (terms_pass2.size() == 1)
1154  return terms_pass2[0].orig;
1155 
1156  // Chop the symmetrized terms into subterms
1157  std::vector<symminfo> sy;
1158  for (auto & i : terms_pass2) {
1159  if (is_exactly_a<add>(i.symm)) {
1160  size_t num = i.symm.nops();
1161  for (size_t j=0; j<num; j++)
1162  sy.push_back(symminfo(i.symm.op(j), i.orig, num));
1163  } else
1164  sy.push_back(symminfo(i.symm, i.orig, 1));
1165  }
1166 
1167  // Sort by symmetrized subterms
1168  std::sort(sy.begin(), sy.end(), symminfo_is_less_by_symmterm());
1169 
1170  // Combine equal symmetrized subterms
1171  std::vector<symminfo> sy_pass2;
1172  exvector result;
1173  for (auto i=sy.begin(); i!=sy.end(); ) {
1174 
1175  // Combine equal terms
1176  auto j = i + 1;
1177  if (j != sy.end() && j->symmterm == i->symmterm) {
1178 
1179  // More than one term, collect the coefficients
1180  ex coeff = i->coeff;
1181  while (j != sy.end() && j->symmterm == i->symmterm) {
1182  coeff += j->coeff;
1183  j++;
1184  }
1185 
1186  // Add combined term to result
1187  if (!coeff.is_zero())
1188  result.push_back(coeff * i->symmterm);
1189 
1190  } else {
1191 
1192  // Single term, store for second pass
1193  sy_pass2.push_back(*i);
1194  }
1195 
1196  i = j;
1197  }
1198 
1199  // Were there any remaining terms that didn't get combined?
1200  if (sy_pass2.size() > 0) {
1201 
1202  // Yes, sort by their original terms
1203  std::sort(sy_pass2.begin(), sy_pass2.end(), symminfo_is_less_by_orig());
1204 
1205  for (std::vector<symminfo>::const_iterator i=sy_pass2.begin(); i!=sy_pass2.end(); ) {
1206 
1207  // How many symmetrized terms of this original term are left?
1208  size_t num = 1;
1209  auto j = i + 1;
1210  while (j != sy_pass2.end() && j->orig == i->orig) {
1211  num++;
1212  j++;
1213  }
1214 
1215  if (num == i->num) {
1216 
1217  // All terms left, then add the original term to the result
1218  result.push_back(i->orig);
1219 
1220  } else {
1221 
1222  // Some terms were combined with others, add up the remaining symmetrized terms
1223  std::vector<symminfo>::const_iterator k;
1224  for (k=i; k!=j; k++)
1225  result.push_back(k->coeff * k->symmterm);
1226  }
1227 
1228  i = j;
1229  }
1230  }
1231 
1232  // Add all resulting terms
1233  ex sum_symm = dynallocate<add>(result);
1234  if (sum_symm.is_zero())
1235  free_indices.clear();
1236  return sum_symm;
1237  }
1238 
1239  // Simplification of products
1240  if (is_exactly_a<mul>(e_expanded)
1241  || is_exactly_a<ncmul>(e_expanded)
1242  || (is_exactly_a<power>(e_expanded) && is_a<indexed>(e_expanded.op(0)) && e_expanded.op(1).is_equal(_ex2)))
1243  return simplify_indexed_product(e_expanded, free_indices, dummy_indices, sp);
1244 
1245  // Cannot do anything
1246  free_indices.clear();
1247  return e_expanded;
1248 }
1249 
1257 {
1258  exvector free_indices, dummy_indices;
1259  scalar_products sp;
1260  return GiNaC::simplify_indexed(*this, free_indices, dummy_indices, sp);
1261 }
1262 
1271 ex ex::simplify_indexed(const scalar_products & sp, unsigned options) const
1272 {
1273  exvector free_indices, dummy_indices;
1274  return GiNaC::simplify_indexed(*this, free_indices, dummy_indices, sp);
1275 }
1276 
1279 {
1280  return GiNaC::symmetrize(*this, get_free_indices());
1281 }
1282 
1285 {
1286  return GiNaC::antisymmetrize(*this, get_free_indices());
1287 }
1288 
1291 {
1292  return GiNaC::symmetrize_cyclic(*this, get_free_indices());
1293 }
1294 
1296 // helper classes
1298 
1299 spmapkey::spmapkey(const ex & v1_, const ex & v2_, const ex & dim_) : dim(dim_)
1300 {
1301  // If indexed, extract base objects
1302  ex s1 = is_a<indexed>(v1_) ? v1_.op(0) : v1_;
1303  ex s2 = is_a<indexed>(v2_) ? v2_.op(0) : v2_;
1304 
1305  // Enforce canonical order in pair
1306  if (s1.compare(s2) > 0) {
1307  v1 = s2;
1308  v2 = s1;
1309  } else {
1310  v1 = s1;
1311  v2 = s2;
1312  }
1313 }
1314 
1315 bool spmapkey::operator==(const spmapkey &other) const
1316 {
1317  if (!v1.is_equal(other.v1))
1318  return false;
1319  if (!v2.is_equal(other.v2))
1320  return false;
1321  if (is_a<wildcard>(dim) || is_a<wildcard>(other.dim))
1322  return true;
1323  else
1324  return dim.is_equal(other.dim);
1325 }
1326 
1327 bool spmapkey::operator<(const spmapkey &other) const
1328 {
1329  int cmp = v1.compare(other.v1);
1330  if (cmp)
1331  return cmp < 0;
1332  cmp = v2.compare(other.v2);
1333  if (cmp)
1334  return cmp < 0;
1335 
1336  // Objects are equal, now check dimensions
1337  if (is_a<wildcard>(dim) || is_a<wildcard>(other.dim))
1338  return false;
1339  else
1340  return dim.compare(other.dim) < 0;
1341 }
1342 
1344 {
1345  std::cerr << "(" << v1 << "," << v2 << "," << dim << ")";
1346 }
1347 
1348 void scalar_products::add(const ex & v1, const ex & v2, const ex & sp)
1349 {
1350  spm[spmapkey(v1, v2)] = sp;
1351 }
1352 
1353 void scalar_products::add(const ex & v1, const ex & v2, const ex & dim, const ex & sp)
1354 {
1355  spm[spmapkey(v1, v2, dim)] = sp;
1356 }
1357 
1358 void scalar_products::add_vectors(const lst & l, const ex & dim)
1359 {
1360  // Add all possible pairs of products
1361  for (auto & it1 : l)
1362  for (auto & it2 : l)
1363  add(it1, it2, it1 * it2);
1364 }
1365 
1367 {
1368  spm.clear();
1369 }
1370 
1372 bool scalar_products::is_defined(const ex & v1, const ex & v2, const ex & dim) const
1373 {
1374  return spm.find(spmapkey(v1, v2, dim)) != spm.end();
1375 }
1376 
1378 ex scalar_products::evaluate(const ex & v1, const ex & v2, const ex & dim) const
1379 {
1380  return spm.find(spmapkey(v1, v2, dim))->second;
1381 }
1382 
1384 {
1385  std::cerr << "map size=" << spm.size() << std::endl;
1386  for (auto & it : spm) {
1387  const spmapkey & k = it.first;
1388  std::cerr << "item key=";
1389  k.debugprint();
1390  std::cerr << ", value=" << it.second << std::endl;
1391  }
1392 }
1393 
1395 {
1396  if (is_a<indexed>(e))
1397  return ex_to<indexed>(e).get_dummy_indices();
1398  else if (is_a<power>(e) && e.op(1)==2) {
1399  return e.op(0).get_free_indices();
1400  }
1401  else if (is_a<mul>(e) || is_a<ncmul>(e)) {
1402  exvector dummies;
1403  exvector free_indices;
1404  for (std::size_t i = 0; i < e.nops(); ++i) {
1405  exvector dummies_of_factor = get_all_dummy_indices_safely(e.op(i));
1406  dummies.insert(dummies.end(), dummies_of_factor.begin(),
1407  dummies_of_factor.end());
1408  exvector free_of_factor = e.op(i).get_free_indices();
1409  free_indices.insert(free_indices.begin(), free_of_factor.begin(),
1410  free_of_factor.end());
1411  }
1412  exvector free_out, dummy_out;
1413  find_free_and_dummy(free_indices.begin(), free_indices.end(), free_out,
1414  dummy_out);
1415  dummies.insert(dummies.end(), dummy_out.begin(), dummy_out.end());
1416  return dummies;
1417  }
1418  else if(is_a<add>(e)) {
1419  exvector result;
1420  for(std::size_t i = 0; i < e.nops(); ++i) {
1421  exvector dummies_of_term = get_all_dummy_indices_safely(e.op(i));
1422  sort(dummies_of_term.begin(), dummies_of_term.end());
1423  exvector new_vec;
1424  set_union(result.begin(), result.end(), dummies_of_term.begin(),
1425  dummies_of_term.end(), std::back_inserter<exvector>(new_vec),
1426  ex_is_less());
1427  result.swap(new_vec);
1428  }
1429  return result;
1430  }
1431  return exvector();
1432 }
1433 
1436 {
1437  exvector p;
1438  bool nc;
1439  product_to_exvector(e, p, nc);
1440  auto ip = p.begin(), ipend = p.end();
1441  exvector v, v1;
1442  while (ip != ipend) {
1443  if (is_a<indexed>(*ip)) {
1444  v1 = ex_to<indexed>(*ip).get_dummy_indices();
1445  v.insert(v.end(), v1.begin(), v1.end());
1446  auto ip1 = ip + 1;
1447  while (ip1 != ipend) {
1448  if (is_a<indexed>(*ip1)) {
1449  v1 = ex_to<indexed>(*ip).get_dummy_indices(ex_to<indexed>(*ip1));
1450  v.insert(v.end(), v1.begin(), v1.end());
1451  }
1452  ++ip1;
1453  }
1454  }
1455  ++ip;
1456  }
1457  return v;
1458 }
1459 
1461 {
1462  exvector common_indices;
1463  set_intersection(va.begin(), va.end(), vb.begin(), vb.end(), std::back_insert_iterator<exvector>(common_indices), ex_is_less());
1464  if (common_indices.empty()) {
1465  return lst{lst{}, lst{}};
1466  } else {
1467  exvector new_indices, old_indices;
1468  old_indices.reserve(2*common_indices.size());
1469  new_indices.reserve(2*common_indices.size());
1470  exvector::const_iterator ip = common_indices.begin(), ipend = common_indices.end();
1471  while (ip != ipend) {
1472  ex newsym = dynallocate<symbol>();
1473  ex newidx;
1474  if(is_exactly_a<spinidx>(*ip))
1475  newidx = dynallocate<spinidx>(newsym, ex_to<spinidx>(*ip).get_dim(),
1476  ex_to<spinidx>(*ip).is_covariant(),
1477  ex_to<spinidx>(*ip).is_dotted());
1478  else if (is_exactly_a<varidx>(*ip))
1479  newidx = dynallocate<varidx>(newsym, ex_to<varidx>(*ip).get_dim(),
1480  ex_to<varidx>(*ip).is_covariant());
1481  else
1482  newidx = dynallocate<idx>(newsym, ex_to<idx>(*ip).get_dim());
1483  old_indices.push_back(*ip);
1484  new_indices.push_back(newidx);
1485  if(is_a<varidx>(*ip)) {
1486  old_indices.push_back(ex_to<varidx>(*ip).toggle_variance());
1487  new_indices.push_back(ex_to<varidx>(newidx).toggle_variance());
1488  }
1489  ++ip;
1490  }
1491  return lst{lst(old_indices.begin(), old_indices.end()), lst(new_indices.begin(), new_indices.end())};
1492  }
1493 }
1494 
1495 ex rename_dummy_indices_uniquely(const exvector & va, const exvector & vb, const ex & b)
1496 {
1497  lst indices_subs = rename_dummy_indices_uniquely(va, vb);
1498  return (indices_subs.op(0).nops()>0 ? b.subs(ex_to<lst>(indices_subs.op(0)), ex_to<lst>(indices_subs.op(1)), subs_options::no_pattern|subs_options::no_index_renaming) : b);
1499 }
1500 
1501 ex rename_dummy_indices_uniquely(const ex & a, const ex & b)
1502 {
1504  if (va.size() > 0) {
1506  if (vb.size() > 0) {
1507  sort(va.begin(), va.end(), ex_is_less());
1508  sort(vb.begin(), vb.end(), ex_is_less());
1509  lst indices_subs = rename_dummy_indices_uniquely(va, vb);
1510  if (indices_subs.op(0).nops() > 0)
1511  return b.subs(ex_to<lst>(indices_subs.op(0)), ex_to<lst>(indices_subs.op(1)), subs_options::no_pattern|subs_options::no_index_renaming);
1512  }
1513  }
1514  return b;
1515 }
1516 
1517 ex rename_dummy_indices_uniquely(exvector & va, const ex & b, bool modify_va)
1518 {
1519  if (va.size() > 0) {
1521  if (vb.size() > 0) {
1522  sort(vb.begin(), vb.end(), ex_is_less());
1523  lst indices_subs = rename_dummy_indices_uniquely(va, vb);
1524  if (indices_subs.op(0).nops() > 0) {
1525  if (modify_va) {
1526  for (auto & i : ex_to<lst>(indices_subs.op(1)))
1527  va.push_back(i);
1528  exvector uncommon_indices;
1529  set_difference(vb.begin(), vb.end(), indices_subs.op(0).begin(), indices_subs.op(0).end(), std::back_insert_iterator<exvector>(uncommon_indices), ex_is_less());
1530  for (auto & ip : uncommon_indices)
1531  va.push_back(ip);
1532  sort(va.begin(), va.end(), ex_is_less());
1533  }
1534  return b.subs(ex_to<lst>(indices_subs.op(0)), ex_to<lst>(indices_subs.op(1)), subs_options::no_pattern|subs_options::no_index_renaming);
1535  }
1536  }
1537  }
1538  return b;
1539 }
1540 
1541 ex expand_dummy_sum(const ex & e, bool subs_idx)
1542 {
1543  ex e_expanded = e.expand();
1545  if (is_a<add>(e_expanded) || is_a<lst>(e_expanded) || is_a<matrix>(e_expanded)) {
1546  return e_expanded.map(fcn);
1547  } else if (is_a<ncmul>(e_expanded) || is_a<mul>(e_expanded) || is_a<power>(e_expanded) || is_a<indexed>(e_expanded)) {
1548  exvector v;
1549  if (is_a<indexed>(e_expanded))
1550  v = ex_to<indexed>(e_expanded).get_dummy_indices();
1551  else
1552  v = get_all_dummy_indices(e_expanded);
1553  ex result = e_expanded;
1554  for (const auto & nu : v) {
1555  if (ex_to<idx>(nu).get_dim().info(info_flags::nonnegint)) {
1556  int idim = ex_to<numeric>(ex_to<idx>(nu).get_dim()).to_int();
1557  ex en = 0;
1558  for (int i=0; i < idim; i++) {
1559  if (subs_idx && is_a<varidx>(nu)) {
1560  ex other = ex_to<varidx>(nu).toggle_variance();
1561  en += result.subs(lst{
1562  nu == idx(i, idim),
1563  other == idx(i, idim)
1564  });
1565  } else {
1566  en += result.subs( nu.op(0) == i );
1567  }
1568  }
1569  result = en;
1570  }
1571  }
1572  return result;
1573  } else {
1574  return e;
1575  }
1576 }
1577 
1578 } // namespace GiNaC
ex idx_symmetrization(const ex &r, const exvector &local_dummy_indices)
Definition: indexed.cpp:762
exvector get_indices() const
Return a vector containing the object&#39;s indices.
Definition: indexed.cpp:423
ex coeff
coefficient of symmetrized term
Definition: indexed.cpp:1011
exvector get_all_dummy_indices(const ex &e)
Returns all dummy indices from the exvector.
Definition: indexed.cpp:1435
bool operator()(const symminfo &si1, const symminfo &si2) const
Definition: indexed.cpp:1018
Interface to GiNaC&#39;s symbolic exponentiation (basis^exponent).
Non-commutative product of expressions.
Definition: ncmul.h:32
Interface to GiNaC&#39;s indexed expressions.
Helper class for storing information about known scalar products which are to be automatically replac...
Definition: indexed.h:227
size_t num
how many symmetrized terms resulted from the original term
Definition: indexed.cpp:1013
int canonicalize(exvector::iterator v, const symmetry &symm)
Canonicalize the order of elements of an expression vector, according to the symmetry properties defi...
Definition: symmetry.cpp:438
This class describes the symmetry of a group of indices.
Definition: symmetry.h:38
ex coeff(const ex &s, int n=1) const
Definition: ex.h:175
idx
Definition: idx.cpp:44
Interface to GiNaC&#39;s symbolic objects.
unsigned hashvalue
hash value
Definition: basic.h:303
ex op(size_t i) const override
Return operand/member at position i.
Definition: container.h:295
exvector get_free_indices() const override
Return a vector containing the free indices of an expression.
Definition: indexed.cpp:494
exvector get_all_dummy_indices_safely(const ex &e)
More reliable version of the form.
Definition: indexed.cpp:1394
void debugprint() const
Definition: indexed.cpp:1343
void print_indexed(const print_context &c, const char *openbrace, const char *closebrace, unsigned level) const
Definition: indexed.cpp:207
int to_int(const numeric &x)
Definition: numeric.h:302
int compare(const ex &other) const
Definition: ex.h:322
Interface to GiNaC&#39;s symmetry definitions.
expands (a+b).i to a.i+b.i
Definition: flags.h:32
ex subs(const exmap &m, unsigned options=0) const
Definition: ex.h:826
void validate() const
Check whether all indices are of class idx and validate the symmetry tree.
Definition: indexed.cpp:368
bool operator()(const symminfo &si1, const symminfo &si2) const
Definition: indexed.cpp:1026
unsigned precedence() const override
Return relative operator precedence (for parenthezing output).
Definition: indexed.h:146
void archive(archive_node &n) const override
Save (a.k.a.
Definition: indexed.cpp:155
bool is_equal(const ex &other) const
Definition: ex.h:345
friend class ex
Definition: basic.h:108
ex expand(unsigned options=0) const
Definition: ex.cpp:73
ex symmetrize_cyclic(const ex &thisex)
Definition: ex.h:805
ex antisymmetrize(const ex &thisex)
Definition: ex.h:799
Definition: add.cpp:38
Product of expressions.
Definition: mul.h:31
ex minimal_dim(const ex &dim1, const ex &dim2)
Return the minimum of two index dimensions.
Definition: idx.cpp:561
This structure stores the original and symmetrized versions of terms obtained during the simplificati...
Definition: indexed.cpp:977
unsigned return_type() const override
Definition: indexed.cpp:323
bool reposition_dummy_indices(ex &e, exvector &variant_dummy_indices, exvector &moved_indices)
Raise/lower dummy indices in a single indexed objects to canonicalize their variance.
Definition: indexed.cpp:618
ex evaluate(const ex &v1, const ex &v2, const ex &dim) const
Return value of defined scalar product pair.
Definition: indexed.cpp:1378
Definition: ex.h:972
Archiving of GiNaC expressions.
Interface to GiNaC&#39;s sums of expressions.
This class is the ABC (abstract base class) of GiNaC&#39;s class hierarchy.
Definition: basic.h:104
size_t nops() const
Definition: ex.h:135
const ex _ex0
Definition: utils.cpp:177
bool are_ex_trivially_equal(const ex &e1, const ex &e2)
Compare two objects of class quickly without doing a deep tree traversal.
Definition: ex.h:684
ex x
Definition: factor.cpp:1641
void add(const ex &v1, const ex &v2, const ex &sp)
Register scalar product pair and its value.
Definition: indexed.cpp:1348
Interface to GiNaC&#39;s products of expressions.
exvector get_free_indices() const override
Return a vector containing the free indices of an expression.
Definition: indexed.cpp:457
ex symm
symmetrized term
Definition: indexed.cpp:982
bool all_index_values_are(unsigned inf) const
Check whether all index values have a certain property.
Definition: indexed.cpp:246
void clear()
Clear all registered scalar products.
Definition: indexed.cpp:1366
Interface to GiNaC&#39;s indices.
container< std::list > lst
Definition: lst.h:32
size_t r
Definition: factor.cpp:770
ex symmetrize_cyclic() const
Symmetrize expression by cyclic permutation over its free indices.
Definition: indexed.cpp:1290
disable pattern matching
Definition: flags.h:51
ex op(size_t i) const
Definition: ex.h:136
Interface to several small and furry utilities needed within GiNaC but not of any interest to the use...
exvector get_dummy_indices() const
Return a vector containing the dummy indices of the object, if any.
Definition: indexed.cpp:429
static void find_variant_indices(const exvector &v, exvector &variant_indices)
Given a set of indices, extract those of class varidx.
Definition: indexed.cpp:602
const ex _ex1
Definition: utils.cpp:193
exvector get_free_indices() const override
Return a vector containing the free indices of an expression.
Definition: indexed.cpp:516
const_iterator begin() const noexcept
Definition: ex.h:647
void read_archive(const archive_node &n, lst &syms) override
Read (a.k.a.
Definition: indexed.cpp:132
Context for latex-parsable output.
Definition: print.h:122
void do_print_tree(const print_tree &c, unsigned level) const
Definition: indexed.cpp:229
size_t number_of_type(const exvector &v)
Definition: indexed.cpp:523
print_func< print_context >(&varidx::do_print). print_func< print_latex >(&varidx
Definition: idx.cpp:45
vector< int > k
Definition: factor.cpp:1466
This class holds an indexed expression.
Definition: indexed.h:39
void find_free_and_dummy(exvector::const_iterator it, exvector::const_iterator itend, exvector &out_free, exvector &out_dummy)
Given a vector of indices, split them into two vectors, one containing the free indices, the other containing the dummy indices (numeric indices are neither free nor dummy ones).
Definition: idx.cpp:521
ex simplify_indexed(const ex &thisex, unsigned options=0)
Definition: ex.h:787
bool has_dummy_index_for(const ex &i) const
Check whether the object has an index that forms a dummy index pair with a given index.
Definition: indexed.cpp:446
terminfo(const ex &orig_, const ex &symm_)
Definition: indexed.cpp:979
exvector get_free_indices() const override
Return a vector containing the free indices of an expression.
Definition: indexed.cpp:479
unsigned options
Definition: factor.cpp:2480
size_t nops() const override
Number of operands/members.
Definition: container.h:118
mvec m
Definition: factor.cpp:771
#define GINAC_ASSERT(X)
Assertion macro for checking invariances.
Definition: assertion.h:33
bool operator()(const ex &lh, const ex &rh) const
Definition: indexed.cpp:723
ex orig
original term
Definition: indexed.cpp:981
ex antisymmetrize() const
Antisymmetrize expression over its free indices.
Definition: indexed.cpp:1284
void do_print(const print_context &c, unsigned level) const
Definition: indexed.cpp:219
GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(add, expairseq, print_func< print_context >(&add::do_print). print_func< print_latex >(&add::do_print_latex). print_func< print_csrc >(&add::do_print_csrc). print_func< print_tree >(&add::do_print_tree). print_func< print_python_repr >(&add::do_print_python_repr)) add
Definition: add.cpp:40
bool is_dummy_pair(const idx &i1, const idx &i2)
Check whether two indices form a dummy pair.
Definition: idx.cpp:502
ex map(map_function &f) const
Definition: ex.h:162
void find_dummy_indices(const exvector &v, exvector &out_dummy)
Given a vector of indices, find the dummy indices.
Definition: idx.h:248
Interface to symbolic matrices.
bool operator()(const ex &e)
Definition: indexed.cpp:510
ex imag_part() const override
Definition: indexed.cpp:306
ex expand_dummy_sum(const ex &e, bool subs_idx)
This function returns the given expression with expanded sums for all dummy index summations...
Definition: indexed.cpp:1541
std::vector< ex > exvector
Definition: basic.h:46
std::map< ex, ex, ex_is_less > exmap
Definition: basic.h:50
void debugprint() const
Definition: indexed.cpp:1383
ex simplify_indexed_product(const ex &e, exvector &free_indices, exvector &dummy_indices, const scalar_products &sp)
Simplify product of indexed expressions (commutative, noncommutative and simple squares), return list of free indices.
Definition: indexed.cpp:779
Interface to GiNaC&#39;s overloaded operators.
ex eval() const override
Perform automatic non-interruptive term rewriting rules.
Definition: indexed.cpp:263
void printindices(const print_context &c, unsigned level) const
Definition: indexed.cpp:165
ex thiscontainer(const exvector &v) const override
Definition: indexed.cpp:313
void shaker_sort(It first, It last, Cmp comp, Swap swapit)
Definition: utils.h:193
void do_print_latex(const print_latex &c, unsigned level) const
Definition: indexed.cpp:224
bool info(unsigned inf) const override
Information about the object.
Definition: indexed.cpp:239
size_t n
Definition: factor.cpp:1463
ex real_part() const override
Definition: indexed.cpp:299
Interface to GiNaC&#39;s symbolic integral.
Base class for print_contexts.
Definition: print.h:102
Definition of GiNaC&#39;s lst.
static ex symm(const ex &e, exvector::const_iterator first, exvector::const_iterator last, bool asymmetric)
Definition: symmetry.cpp:485
bool operator()(const terminfo &ti1, const terminfo &ti2) const
Definition: indexed.cpp:987
ex op(size_t i) const override
Return operand/member at position i.
This structure stores the individual symmetrized terms obtained during the simplification of sums...
Definition: indexed.cpp:995
static void product_to_exvector(const ex &e, exvector &v, bool &non_commutative)
Definition: indexed.cpp:731
indexed(const ex &b)
Construct indexed object with no index.
Definition: indexed.cpp:64
ex derivative(const symbol &s) const override
Implementation of ex::diff() for an indexed object always returns 0.
Definition: indexed.cpp:388
bool operator()(const ex &lh, const ex &rh) const
Definition: indexed.cpp:398
ex simplify_indexed(unsigned options=0) const
Simplify/canonicalize expression containing indexed objects.
Definition: indexed.cpp:1256
virtual ex eval_indexed(const basic &i) const
Perform automatic symbolic evaluations on indexed expression that contains this object as the base ex...
Definition: basic.cpp:465
This class stores all properties needed to record/retrieve the state of one object of class basic (or...
Definition: archive.h:48
Lightweight wrapper for GiNaC&#39;s symbolic objects.
Definition: ex.h:72
symminfo(const ex &symmterm_, const ex &orig_, size_t num_)
Definition: indexed.cpp:999
container< std::vector > exprseq
Definition: exprseq.h:32
ex symmetrize() const
Symmetrize expression over its free indices.
Definition: indexed.cpp:1278
symmetry sy_symm()
Definition: symmetry.h:121
Interface to GiNaC&#39;s non-commutative products of expressions.
lst rename_dummy_indices_uniquely(const exvector &va, const exvector &vb)
Similar to above, where va and vb are the same and the return value is a list of two lists for substi...
Definition: indexed.cpp:1460
ex symmetrize(const ex &thisex)
Definition: ex.h:793
Interface to GiNaC&#39;s initially known functions.
ex coeff(const ex &thisex, const ex &s, int n=1)
Definition: ex.h:742
void validate(unsigned n)
Verify that all indices of this node are in the range [0..n-1].
Definition: symmetry.cpp:312
Basic CAS symbol.
Definition: symbol.h:38
virtual int compare_same_type(const basic &other) const
Returns order relation between two objects of same type.
Definition: basic.cpp:719
unsigned return_type() const
Definition: ex.h:230
exvector get_free_indices() const override
Return a vector containing the free indices of an expression.
Definition: indexed.cpp:464
Wrapper template for making GiNaC classes out of STL containers.
Definition: container.h:73
bool operator==(const spmapkey &other) const
Definition: indexed.cpp:1315
bool is_defined(const ex &v1, const ex &v2, const ex &dim) const
Check whether scalar product pair is defined.
Definition: indexed.cpp:1372
Interface to relations between expressions.
static ex rename_dummy_indices(const ex &e, exvector &global_dummy_indices, exvector &local_dummy_indices)
Rename dummy indices in an expression.
Definition: indexed.cpp:540
bool is_zero() const
Definition: ex.h:213
ex symmterm
symmetrized term
Definition: indexed.cpp:1010
void add_vectors(const lst &l, const ex &dim=wild())
Register list of vectors.
Definition: indexed.cpp:1358
symmetry sy_anti()
Definition: symmetry.h:126
void reserve(size_t)
Definition: container.h:54
static bool indices_consistent(const exvector &v1, const exvector &v2)
Check whether two sorted index vectors are consistent (i.e.
Definition: indexed.cpp:414
size_t c
Definition: factor.cpp:770
const symmetry & not_symmetric()
Definition: symmetry.cpp:350
size_t nops() const override
Number of operands/members.
ex orig
original term
Definition: indexed.cpp:1012
ex expand(unsigned options=0) const override
Expand expression, i.e.
Definition: indexed.cpp:331
bool operator<(const spmapkey &other) const
Definition: indexed.cpp:1327
Context for tree-like output for debugging.
Definition: print.h:146
ex symtree
Index symmetry (tree of symmetry objects)
Definition: indexed.h:202
unsigned flags
of type status_flags
Definition: basic.h:302
bool hasindex(const ex &x, const ex &sym)
Definition: indexed.cpp:1032
ex expand(const ex &thisex, unsigned options=0)
Definition: ex.h:715
const ex _ex2
Definition: utils.cpp:197
GINAC_BIND_UNARCHIVER(add)
exvector get_free_indices() const
Definition: ex.h:206

This page is part of the GiNaC developer's reference. It was generated automatically by doxygen. For an introduction, see the tutorial.