GiNaC  1.8.3
indexed.cpp
Go to the documentation of this file.
1 
5 /*
6  * GiNaC Copyright (C) 1999-2022 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
Interface to GiNaC's sums of expressions.
Archiving of GiNaC expressions.
#define GINAC_ASSERT(X)
Assertion macro for checking invariances.
Definition: assertion.h:33
exvector get_free_indices() const override
Return a vector containing the free indices of an expression.
Definition: indexed.cpp:464
This class stores all properties needed to record/retrieve the state of one object of class basic (or...
Definition: archive.h:49
This class is the ABC (abstract base class) of GiNaC's class hierarchy.
Definition: basic.h:105
unsigned hashvalue
hash value
Definition: basic.h:303
friend class ex
Definition: basic.h:108
unsigned flags
of type status_flags
Definition: basic.h:302
virtual int compare_same_type(const basic &other) const
Returns order relation between two objects of same type.
Definition: basic.cpp:719
void reserve(size_t)
Definition: container.h:54
Wrapper template for making GiNaC classes out of STL containers.
Definition: container.h:73
size_t nops() const override
Number of operands/members.
Definition: container.h:118
ex op(size_t i) const override
Return operand/member at position i.
Definition: container.h:295
Lightweight wrapper for GiNaC's symbolic objects.
Definition: ex.h:72
ex map(map_function &f) const
Definition: ex.h:162
const_iterator begin() const noexcept
Definition: ex.h:647
exvector get_free_indices() const
Definition: ex.h:206
ex expand(unsigned options=0) const
Definition: ex.cpp:73
bool is_equal(const ex &other) const
Definition: ex.h:345
ex simplify_indexed(unsigned options=0) const
Simplify/canonicalize expression containing indexed objects.
Definition: indexed.cpp:1256
ex symmetrize_cyclic() const
Symmetrize expression by cyclic permutation over its free indices.
Definition: indexed.cpp:1290
unsigned return_type() const
Definition: ex.h:230
const_iterator end() const noexcept
Definition: ex.h:652
size_t nops() const
Definition: ex.h:135
ex subs(const exmap &m, unsigned options=0) const
Definition: ex.h:826
int compare(const ex &other) const
Definition: ex.h:322
ex symmetrize() const
Symmetrize expression over its free indices.
Definition: indexed.cpp:1278
bool is_zero() const
Definition: ex.h:213
ex op(size_t i) const
Definition: ex.h:136
ex coeff(const ex &s, int n=1) const
Definition: ex.h:175
ex antisymmetrize() const
Antisymmetrize expression over its free indices.
Definition: indexed.cpp:1284
size_t nops() const override
Number of operands/members.
ex op(size_t i) const override
Return operand/member at position i.
@ expand_indexed
expands (a+b).i to a.i+b.i
Definition: flags.h:32
This class holds an indexed expression.
Definition: indexed.h:40
bool all_index_values_are(unsigned inf) const
Check whether all index values have a certain property.
Definition: indexed.cpp:246
ex thiscontainer(const exvector &v) const override
Definition: indexed.cpp:313
unsigned precedence() const override
Return relative operator precedence (for parenthezing output).
Definition: indexed.h:146
void printindices(const print_context &c, unsigned level) const
Definition: indexed.cpp:165
exvector get_free_indices() const override
Return a vector containing the free indices of an expression.
Definition: indexed.cpp:457
void archive(archive_node &n) const override
Save (a.k.a.
Definition: indexed.cpp:155
void do_print(const print_context &c, unsigned level) const
Definition: indexed.cpp:219
ex expand(unsigned options=0) const override
Expand expression, i.e.
Definition: indexed.cpp:331
void do_print_latex(const print_latex &c, unsigned level) const
Definition: indexed.cpp:224
void read_archive(const archive_node &n, lst &syms) override
Read (a.k.a.
Definition: indexed.cpp:132
void do_print_tree(const print_tree &c, unsigned level) const
Definition: indexed.cpp:229
ex derivative(const symbol &s) const override
Implementation of ex::diff() for an indexed object always returns 0.
Definition: indexed.cpp:388
bool info(unsigned inf) const override
Information about the object.
Definition: indexed.cpp:239
ex eval() const override
Perform automatic non-interruptive term rewriting rules.
Definition: indexed.cpp:263
indexed(const ex &b)
Construct indexed object with no index.
Definition: indexed.cpp:64
ex real_part() const override
Definition: indexed.cpp:299
void validate() const
Check whether all indices are of class idx and validate the symmetry tree.
Definition: indexed.cpp:368
ex symtree
Index symmetry (tree of symmetry objects)
Definition: indexed.h:202
exvector get_indices() const
Return a vector containing the object's indices.
Definition: indexed.cpp:423
exvector get_dummy_indices() const
Return a vector containing the dummy indices of the object, if any.
Definition: indexed.cpp:429
void print_indexed(const print_context &c, const char *openbrace, const char *closebrace, unsigned level) const
Definition: indexed.cpp:207
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
unsigned return_type() const override
Definition: indexed.cpp:323
ex imag_part() const override
Definition: indexed.cpp:306
exvector get_free_indices() const override
Return a vector containing the free indices of an expression.
Definition: indexed.cpp:516
Product of expressions.
Definition: mul.h:32
exvector get_free_indices() const override
Return a vector containing the free indices of an expression.
Definition: indexed.cpp:479
Non-commutative product of expressions.
Definition: ncmul.h:33
exvector get_free_indices() const override
Return a vector containing the free indices of an expression.
Definition: indexed.cpp:494
Base class for print_contexts.
Definition: print.h:103
Context for latex-parsable output.
Definition: print.h:123
Context for tree-like output for debugging.
Definition: print.h:147
Helper class for storing information about known scalar products which are to be automatically replac...
Definition: indexed.h:227
bool is_defined(const ex &v1, const ex &v2, const ex &dim) const
Check whether scalar product pair is defined.
Definition: indexed.cpp:1372
void debugprint() const
Definition: indexed.cpp:1383
ex evaluate(const ex &v1, const ex &v2, const ex &dim) const
Return value of defined scalar product pair.
Definition: indexed.cpp:1378
void add(const ex &v1, const ex &v2, const ex &sp)
Register scalar product pair and its value.
Definition: indexed.cpp:1348
void clear()
Clear all registered scalar products.
Definition: indexed.cpp:1366
void add_vectors(const lst &l, const ex &dim=wild())
Register list of vectors.
Definition: indexed.cpp:1358
void debugprint() const
Definition: indexed.cpp:1343
bool operator<(const spmapkey &other) const
Definition: indexed.cpp:1327
bool operator==(const spmapkey &other) const
Definition: indexed.cpp:1315
@ no_pattern
disable pattern matching
Definition: flags.h:51
Basic CAS symbol.
Definition: symbol.h:39
This class describes the symmetry of a group of indices.
Definition: symmetry.h:39
bool operator()(const symminfo &si1, const symminfo &si2) const
Definition: indexed.cpp:1026
bool operator()(const symminfo &si1, const symminfo &si2) const
Definition: indexed.cpp:1018
This structure stores the individual symmetrized terms obtained during the simplification of sums.
Definition: indexed.cpp:995
ex coeff
coefficient of symmetrized term
Definition: indexed.cpp:1011
ex symmterm
symmetrized term
Definition: indexed.cpp:1010
ex orig
original term
Definition: indexed.cpp:1012
size_t num
how many symmetrized terms resulted from the original term
Definition: indexed.cpp:1013
symminfo(const ex &symmterm_, const ex &orig_, size_t num_)
Definition: indexed.cpp:999
bool operator()(const terminfo &ti1, const terminfo &ti2) const
Definition: indexed.cpp:987
This structure stores the original and symmetrized versions of terms obtained during the simplificati...
Definition: indexed.cpp:977
terminfo(const ex &orig_, const ex &symm_)
Definition: indexed.cpp:979
ex symm
symmetrized term
Definition: indexed.cpp:982
ex orig
original term
Definition: indexed.cpp:981
unsigned options
Definition: factor.cpp:2480
vector< int > k
Definition: factor.cpp:1466
ex x
Definition: factor.cpp:1641
size_t n
Definition: factor.cpp:1463
size_t c
Definition: factor.cpp:770
size_t r
Definition: factor.cpp:770
mvec m
Definition: factor.cpp:771
Interface to GiNaC's indices.
Interface to GiNaC's indexed expressions.
Interface to GiNaC's initially known functions.
Interface to GiNaC's symbolic integral.
Definition of GiNaC's lst.
Interface to symbolic matrices.
Interface to GiNaC's products of expressions.
Definition: add.cpp:38
static bool indices_consistent(const exvector &v1, const exvector &v2)
Check whether two sorted index vectors are consistent (i.e.
Definition: indexed.cpp:414
ex idx_symmetrization(const ex &r, const exvector &local_dummy_indices)
Definition: indexed.cpp:762
const ex _ex2
Definition: utils.cpp:389
ex minimal_dim(const ex &dim1, const ex &dim2)
Return the minimum of two index dimensions.
Definition: idx.cpp:561
const symmetry & not_symmetric()
Definition: symmetry.cpp:350
container< std::list > lst
Definition: lst.h:32
std::map< ex, ex, ex_is_less > exmap
Definition: basic.h:50
symmetry sy_symm()
Definition: symmetry.h:121
ex symmetrize(const ex &thisex)
Definition: ex.h:793
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:385
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
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
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
ex simplify_indexed(const ex &thisex, unsigned options=0)
Definition: ex.h:787
idx
Definition: idx.cpp:44
bool is_dummy_pair(const idx &i1, const idx &i2)
Check whether two indices form a dummy pair.
Definition: idx.cpp:502
print_func< print_context >(&varidx::do_print). print_func< print_latex >(&varidx
Definition: idx.cpp:45
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
GINAC_IMPLEMENT_REGISTERED_CLASS_OPT_T(lst, basic, print_func< print_context >(&lst::do_print). print_func< print_tree >(&lst::do_print_tree)) template<> bool lst GINAC_BIND_UNARCHIVER(lst)
Specialization of container::info() for lst.
Definition: lst.cpp:42
symmetry sy_anti()
Definition: symmetry.h:126
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),...
Definition: indexed.cpp:779
ex antisymmetrize(const ex &thisex)
Definition: ex.h:799
ex op(const ex &thisex, size_t i)
Definition: ex.h:811
void shaker_sort(It first, It last, Cmp comp, Swap swapit)
Definition: utils.h:193
ex coeff(const ex &thisex, const ex &s, int n=1)
Definition: ex.h:742
static ex symm(const ex &e, exvector::const_iterator first, exvector::const_iterator last, bool asymmetric)
Definition: symmetry.cpp:485
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
static void product_to_exvector(const ex &e, exvector &v, bool &non_commutative)
Definition: indexed.cpp:731
size_t number_of_type(const exvector &v)
Definition: indexed.cpp:523
void find_dummy_indices(const exvector &v, exvector &out_dummy)
Given a vector of indices, find the dummy indices.
Definition: idx.h:248
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_cyclic(const ex &thisex)
Definition: ex.h:805
const ex _ex0
Definition: utils.cpp:369
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,...
Definition: idx.cpp:521
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
bool hasindex(const ex &x, const ex &sym)
Definition: indexed.cpp:1032
exvector get_all_dummy_indices(const ex &e)
Returns all dummy indices from the exvector.
Definition: indexed.cpp:1435
container< std::vector > exprseq
Definition: exprseq.h:32
exvector get_all_dummy_indices_safely(const ex &e)
More reliable version of the form.
Definition: indexed.cpp:1394
ex expand(const ex &thisex, unsigned options=0)
Definition: ex.h:715
Definition: ex.h:972
Interface to GiNaC's non-commutative products of expressions.
Interface to GiNaC's overloaded operators.
Interface to GiNaC's symbolic exponentiation (basis^exponent).
Interface to relations between expressions.
bool operator()(const ex &lh, const ex &rh) const
Definition: indexed.cpp:723
bool operator()(const ex &lh, const ex &rh) const
Definition: indexed.cpp:398
bool operator()(const ex &e)
Definition: indexed.cpp:510
Interface to GiNaC's symbolic objects.
Interface to GiNaC's symmetry definitions.
Interface to several small and furry utilities needed within GiNaC but not of any interest to the use...

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