ginac/tensor.cpp

   1 /** @file tensor.cpp
   2  *
   3  *  Implementation of GiNaC's special tensors. */
   4
   5 /*
   6  *  GiNaC Copyright (C) 1999-2008 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 <iostream>
  24 #include <stdexcept>
  25 #include <vector>
  26
  27 #include "tensor.h"
  28 #include "idx.h"
  29 #include "indexed.h"
  30 #include "symmetry.h"
  31 #include "relational.h"
  32 #include "operators.h"
  33 #include "lst.h"
  34 #include "numeric.h"
  35 #include "matrix.h"
  36 #include "archive.h"
  37 #include "utils.h"
  38
  39 namespace GiNaC {
  40
  41 GINAC_IMPLEMENT_REGISTERED_CLASS(tensor, basic)
  42
  43 GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(tensdelta, tensor,
  44   print_func<print_dflt>(&tensdelta::do_print).
  45   print_func<print_latex>(&tensdelta::do_print_latex))
  46
  47 GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(tensmetric, tensor,
  48   print_func<print_dflt>(&tensmetric::do_print).
  49   print_func<print_latex>(&tensmetric::do_print))
  50
  51 GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(minkmetric, tensmetric,
  52   print_func<print_dflt>(&minkmetric::do_print).
  53   print_func<print_latex>(&minkmetric::do_print_latex))
  54
  55 GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(spinmetric, tensmetric,
  56   print_func<print_dflt>(&spinmetric::do_print).
  57   print_func<print_latex>(&spinmetric::do_print_latex))
  58
  59 GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(tensepsilon, tensor,
  60   print_func<print_dflt>(&tensepsilon::do_print).
  61   print_func<print_latex>(&tensepsilon::do_print_latex))
  62
  63 //////////
  64 // constructors
  65 //////////
  66
  67 tensor::tensor()
  68 {
  69         setflag(status_flags::evaluated | status_flags::expanded);
  70 }
  71
  72 DEFAULT_CTOR(tensdelta)
  73 DEFAULT_CTOR(tensmetric)
  74
  75 minkmetric::minkmetric() : pos_sig(false)
  76 {
  77 }
  78
  79 spinmetric::spinmetric()
  80 {
  81 }
  82
  83 minkmetric::minkmetric(bool ps) : pos_sig(ps)
  84 {
  85 }
  86
  87 tensepsilon::tensepsilon() : minkowski(false), pos_sig(false)
  88 {
  89 }
  90
  91 tensepsilon::tensepsilon(bool mink, bool ps) : minkowski(mink), pos_sig(ps)
  92 {
  93 }
  94
  95 //////////
  96 // archiving
  97 //////////
  98
  99 DEFAULT_ARCHIVING(tensor)
 100 DEFAULT_ARCHIVING(tensdelta)
 101 DEFAULT_ARCHIVING(tensmetric)
 102 DEFAULT_ARCHIVING(spinmetric)
 103 DEFAULT_UNARCHIVE(minkmetric)
 104 DEFAULT_UNARCHIVE(tensepsilon)
 105
 106 minkmetric::minkmetric(const archive_node &n, lst &sym_lst) : inherited(n, sym_lst)
 107 {
 108         n.find_bool("pos_sig", pos_sig);
 109 }
 110
 111 void minkmetric::archive(archive_node &n) const
 112 {
 113         inherited::archive(n);
 114         n.add_bool("pos_sig", pos_sig);
 115 }
 116
 117 tensepsilon::tensepsilon(const archive_node &n, lst &sym_lst) : inherited(n, sym_lst)
 118 {
 119         n.find_bool("minkowski", minkowski);
 120         n.find_bool("pos_sig", pos_sig);
 121 }
 122
 123 void tensepsilon::archive(archive_node &n) const
 124 {
 125         inherited::archive(n);
 126         n.add_bool("minkowski", minkowski);
 127         n.add_bool("pos_sig", pos_sig);
 128 }
 129
 130 //////////
 131 // functions overriding virtual functions from base classes
 132 //////////
 133
 134 DEFAULT_COMPARE(tensor)
 135 DEFAULT_COMPARE(tensdelta)
 136 DEFAULT_COMPARE(tensmetric)
 137 DEFAULT_COMPARE(spinmetric)
 138
 139 bool tensdelta::info(unsigned inf) const
 140 {
 141         if(inf == info_flags::real)
 142                 return true;
 143
 144         return false;
 145 }
 146
 147 bool tensmetric::info(unsigned inf) const
 148 {
 149         if(inf == info_flags::real)
 150                 return true;
 151
 152         return false;
 153 }
 154
 155 int minkmetric::compare_same_type(const basic & other) const
 156 {
 157         GINAC_ASSERT(is_a<minkmetric>(other));
 158         const minkmetric &o = static_cast<const minkmetric &>(other);
 159
 160         if (pos_sig != o.pos_sig)
 161                 return pos_sig ? -1 : 1;
 162         else
 163                 return inherited::compare_same_type(other);
 164 }
 165
 166 bool minkmetric::info(unsigned inf) const
 167 {
 168         if(inf == info_flags::real)
 169                 return true;
 170
 171         return false;
 172 }
 173
 174 int tensepsilon::compare_same_type(const basic & other) const
 175 {
 176         GINAC_ASSERT(is_a<tensepsilon>(other));
 177         const tensepsilon &o = static_cast<const tensepsilon &>(other);
 178
 179         if (minkowski != o.minkowski)
 180                 return minkowski ? -1 : 1;
 181         else if (pos_sig != o.pos_sig)
 182                 return pos_sig ? -1 : 1;
 183         else
 184                 return inherited::compare_same_type(other);
 185 }
 186
 187 bool tensepsilon::info(unsigned inf) const
 188 {
 189         if(inf == info_flags::real)
 190                 return true;
 191
 192         return false;
 193 }
 194
 195 bool spinmetric::info(unsigned inf) const
 196 {
 197         if(inf == info_flags::real)
 198                 return true;
 199
 200         return false;
 201 }
 202
 203 DEFAULT_PRINT_LATEX(tensdelta, "delta", "\\delta")
 204 DEFAULT_PRINT(tensmetric, "g")
 205 DEFAULT_PRINT_LATEX(minkmetric, "eta", "\\eta")
 206 DEFAULT_PRINT_LATEX(spinmetric, "eps", "\\varepsilon")
 207 DEFAULT_PRINT_LATEX(tensepsilon, "eps", "\\varepsilon")
 208
 209 /** Automatic symbolic evaluation of an indexed delta tensor. */
 210 ex tensdelta::eval_indexed(const basic & i) const
 211 {
 212         GINAC_ASSERT(is_a<indexed>(i));
 213         GINAC_ASSERT(i.nops() == 3);
 214         GINAC_ASSERT(is_a<tensdelta>(i.op(0)));
 215
 216         const idx & i1 = ex_to<idx>(i.op(1));
 217         const idx & i2 = ex_to<idx>(i.op(2));
 218
 219         // The dimension of the indices must be equal, otherwise we use the minimal
 220         // dimension
 221         if (!i1.get_dim().is_equal(i2.get_dim())) {
 222                 ex min_dim = i1.minimal_dim(i2);
 223                 exmap m;
 224                 m[i1] = i1.replace_dim(min_dim);
 225                 m[i2] = i2.replace_dim(min_dim);
 226                 return i.subs(m, subs_options::no_pattern);
 227         }
 228
 229         // Trace of delta tensor is the (effective) dimension of the space
 230         if (is_dummy_pair(i1, i2)) {
 231                 try {
 232                         return i1.minimal_dim(i2);
 233                 } catch (std::exception &e) {
 234                         return i.hold();
 235                 }
 236         }
 237
 238         // Numeric evaluation
 239         if (static_cast<const indexed &>(i).all_index_values_are(info_flags::integer)) {
 240                 int n1 = ex_to<numeric>(i1.get_value()).to_int(), n2 = ex_to<numeric>(i2.get_value()).to_int();
 241                 if (n1 == n2)
 242                         return _ex1;
 243                 else
 244                         return _ex0;
 245         }
 246
 247         // No further simplifications
 248         return i.hold();
 249 }
 250
 251 /** Automatic symbolic evaluation of an indexed metric tensor. */
 252 ex tensmetric::eval_indexed(const basic & i) const
 253 {
 254         GINAC_ASSERT(is_a<indexed>(i));
 255         GINAC_ASSERT(i.nops() == 3);
 256         GINAC_ASSERT(is_a<tensmetric>(i.op(0)));
 257         GINAC_ASSERT(is_a<varidx>(i.op(1)));
 258         GINAC_ASSERT(is_a<varidx>(i.op(2)));
 259
 260         const varidx & i1 = ex_to<varidx>(i.op(1));
 261         const varidx & i2 = ex_to<varidx>(i.op(2));
 262
 263         // The dimension of the indices must be equal, otherwise we use the minimal
 264         // dimension
 265         if (!i1.get_dim().is_equal(i2.get_dim())) {
 266                 ex min_dim = i1.minimal_dim(i2);
 267                 exmap m;
 268                 m[i1] = i1.replace_dim(min_dim);
 269                 m[i2] = i2.replace_dim(min_dim);
 270                 return i.subs(m, subs_options::no_pattern);
 271         }
 272
 273         // A metric tensor with one covariant and one contravariant index gets
 274         // replaced by a delta tensor
 275         if (i1.is_covariant() != i2.is_covariant())
 276                 return delta_tensor(i1, i2);
 277
 278         // No further simplifications
 279         return i.hold();
 280 }
 281
 282 /** Automatic symbolic evaluation of an indexed Lorentz metric tensor. */
 283 ex minkmetric::eval_indexed(const basic & i) const
 284 {
 285         GINAC_ASSERT(is_a<indexed>(i));
 286         GINAC_ASSERT(i.nops() == 3);
 287         GINAC_ASSERT(is_a<minkmetric>(i.op(0)));
 288         GINAC_ASSERT(is_a<varidx>(i.op(1)));
 289         GINAC_ASSERT(is_a<varidx>(i.op(2)));
 290
 291         const varidx & i1 = ex_to<varidx>(i.op(1));
 292         const varidx & i2 = ex_to<varidx>(i.op(2));
 293
 294         // Numeric evaluation
 295         if (static_cast<const indexed &>(i).all_index_values_are(info_flags::nonnegint)) {
 296                 int n1 = ex_to<numeric>(i1.get_value()).to_int(), n2 = ex_to<numeric>(i2.get_value()).to_int();
 297                 if (n1 != n2)
 298                         return _ex0;
 299                 else if (n1 == 0)
 300                         return pos_sig ? _ex_1 : _ex1;
 301                 else
 302                         return pos_sig ? _ex1 : _ex_1;
 303         }
 304
 305         // Perform the usual evaluations of a metric tensor
 306         return inherited::eval_indexed(i);
 307 }
 308
 309 /** Automatic symbolic evaluation of an indexed metric tensor. */
 310 ex spinmetric::eval_indexed(const basic & i) const
 311 {
 312         GINAC_ASSERT(is_a<indexed>(i));
 313         GINAC_ASSERT(i.nops() == 3);
 314         GINAC_ASSERT(is_a<spinmetric>(i.op(0)));
 315         GINAC_ASSERT(is_a<spinidx>(i.op(1)));
 316         GINAC_ASSERT(is_a<spinidx>(i.op(2)));
 317
 318         const spinidx & i1 = ex_to<spinidx>(i.op(1));
 319         const spinidx & i2 = ex_to<spinidx>(i.op(2));
 320
 321         // Convolutions are zero
 322         if (!(static_cast<const indexed &>(i).get_dummy_indices().empty()))
 323                 return _ex0;
 324
 325         // Numeric evaluation
 326         if (static_cast<const indexed &>(i).all_index_values_are(info_flags::nonnegint)) {
 327                 int n1 = ex_to<numeric>(i1.get_value()).to_int(), n2 = ex_to<numeric>(i2.get_value()).to_int();
 328                 if (n1 == n2)
 329                         return _ex0;
 330                 else if (n1 < n2)
 331                         return _ex1;
 332                 else
 333                         return _ex_1;
 334         }
 335
 336         // No further simplifications
 337         return i.hold();
 338 }
 339
 340 /** Automatic symbolic evaluation of an indexed epsilon tensor. */
 341 ex tensepsilon::eval_indexed(const basic & i) const
 342 {
 343         GINAC_ASSERT(is_a<indexed>(i));
 344         GINAC_ASSERT(i.nops() > 1);
 345         GINAC_ASSERT(is_a<tensepsilon>(i.op(0)));
 346
 347         // Convolutions are zero
 348         if (!(static_cast<const indexed &>(i).get_dummy_indices().empty()))
 349                 return _ex0;
 350
 351         // Numeric evaluation
 352         if (static_cast<const indexed &>(i).all_index_values_are(info_flags::nonnegint)) {
 353
 354                 // Get sign of index permutation (the indices should already be in
 355                 // a canonic order but we can't assume what exactly that order is)
 356                 std::vector<int> v;
 357                 v.reserve(i.nops() - 1);
 358                 for (size_t j=1; j<i.nops(); j++)
 359                         v.push_back(ex_to<numeric>(ex_to<idx>(i.op(j)).get_value()).to_int());
 360                 int sign = permutation_sign(v.begin(), v.end());
 361
 362                 // In a Minkowski space, check for covariant indices
 363                 if (minkowski) {
 364                         for (size_t j=1; j<i.nops(); j++) {
 365                                 const ex & x = i.op(j);
 366                                 if (!is_a<varidx>(x))
 367                                         throw(std::runtime_error("indices of epsilon tensor in Minkowski space must be of type varidx"));
 368                                 if (ex_to<varidx>(x).is_covariant())
 369                                         if (ex_to<idx>(x).get_value().is_zero())
 370                                                 sign = (pos_sig ? -sign : sign);
 371                                         else
 372                                                 sign = (pos_sig ? sign : -sign);
 373                         }
 374                 }
 375
 376                 return sign;
 377         }
 378
 379         // No further simplifications
 380         return i.hold();
 381 }
 382
 383 bool tensor::replace_contr_index(exvector::iterator self, exvector::iterator other) const
 384 {
 385         // Try to contract the first index
 386         const idx *self_idx = &ex_to<idx>(self->op(1));
 387         const idx *free_idx = &ex_to<idx>(self->op(2));
 388         bool first_index_tried = false;
 389
 390 again:
 391         if (self_idx->is_symbolic()) {
 392                 for (size_t i=1; i<other->nops(); i++) {
 393                         if (! is_a<idx>(other->op(i)))
 394                                 continue;
 395                         const idx &other_idx = ex_to<idx>(other->op(i));
 396                         if (is_dummy_pair(*self_idx, other_idx)) {
 397
 398                                 // Contraction found, remove this tensor and substitute the
 399                                 // index in the second object
 400                                 try {
 401                                         // minimal_dim() throws an exception when index dimensions are not comparable
 402                                         ex min_dim = self_idx->minimal_dim(other_idx);
 403                                         *other = other->subs(other_idx == free_idx->replace_dim(min_dim));
 404                                         *self = _ex1; // *other is assigned first because assigning *self invalidates free_idx
 405                                         return true;
 406                                 } catch (std::exception &e) {
 407                                         return false;
 408                                 }
 409                         }
 410                 }
 411         }
 412
 413         if (!first_index_tried) {
 414
 415                 // No contraction with the first index found, try the second index
 416                 self_idx = &ex_to<idx>(self->op(2));
 417                 free_idx = &ex_to<idx>(self->op(1));
 418                 first_index_tried = true;
 419                 goto again;
 420         }
 421
 422         return false;
 423 }
 424
 425 /** Contraction of an indexed delta tensor with something else. */
 426 bool tensdelta::contract_with(exvector::iterator self, exvector::iterator other, exvector & v) const
 427 {
 428         GINAC_ASSERT(is_a<indexed>(*self));
 429         GINAC_ASSERT(is_a<indexed>(*other));
 430         GINAC_ASSERT(self->nops() == 3);
 431         GINAC_ASSERT(is_a<tensdelta>(self->op(0)));
 432
 433         // Replace the dummy index with this tensor's other index and remove
 434         // the tensor (this is valid for contractions with all other tensors)
 435         return replace_contr_index(self, other);
 436 }
 437
 438 /** Contraction of an indexed metric tensor with something else. */
 439 bool tensmetric::contract_with(exvector::iterator self, exvector::iterator other, exvector & v) const
 440 {
 441         GINAC_ASSERT(is_a<indexed>(*self));
 442         GINAC_ASSERT(is_a<indexed>(*other));
 443         GINAC_ASSERT(self->nops() == 3);
 444         GINAC_ASSERT(is_a<tensmetric>(self->op(0)));
 445
 446         // If contracting with the delta tensor, let the delta do it
 447         // (don't raise/lower delta indices)
 448         if (is_a<tensdelta>(other->op(0)))
 449                 return false;
 450
 451         // Replace the dummy index with this tensor's other index and remove
 452         // the tensor
 453         return replace_contr_index(self, other);
 454 }
 455
 456 /** Contraction of an indexed spinor metric with something else. */
 457 bool spinmetric::contract_with(exvector::iterator self, exvector::iterator other, exvector & v) const
 458 {
 459         GINAC_ASSERT(is_a<indexed>(*self));
 460         GINAC_ASSERT(is_a<indexed>(*other));
 461         GINAC_ASSERT(self->nops() == 3);
 462         GINAC_ASSERT(is_a<spinmetric>(self->op(0)));
 463
 464         // Contractions between spinor metrics
 465         if (is_a<spinmetric>(other->op(0))) {
 466                 const idx &self_i1 = ex_to<idx>(self->op(1));
 467                 const idx &self_i2 = ex_to<idx>(self->op(2));
 468                 const idx &other_i1 = ex_to<idx>(other->op(1));
 469                 const idx &other_i2 = ex_to<idx>(other->op(2));
 470
 471                 if (is_dummy_pair(self_i1, other_i1)) {
 472                         if (is_dummy_pair(self_i2, other_i2))
 473                                 *self = _ex2;
 474                         else
 475                                 *self = delta_tensor(self_i2, other_i2);
 476                         *other = _ex1;
 477                         return true;
 478                 } else if (is_dummy_pair(self_i1, other_i2)) {
 479                         if (is_dummy_pair(self_i2, other_i1))
 480                                 *self = _ex_2;
 481                         else
 482                                 *self = -delta_tensor(self_i2, other_i1);
 483                         *other = _ex1;
 484                         return true;
 485                 } else if (is_dummy_pair(self_i2, other_i1)) {
 486                         *self = -delta_tensor(self_i1, other_i2);
 487                         *other = _ex1;
 488                         return true;
 489                 } else if (is_dummy_pair(self_i2, other_i2)) {
 490                         *self = delta_tensor(self_i1, other_i1);
 491                         *other = _ex1;
 492                         return true;
 493                 }
 494         }
 495
 496         // If contracting with the delta tensor, let the delta do it
 497         // (don't raise/lower delta indices)
 498         if (is_a<tensdelta>(other->op(0)))
 499                 return false;
 500
 501         // Try to contract first index
 502         const idx *self_idx = &ex_to<idx>(self->op(1));
 503         const idx *free_idx = &ex_to<idx>(self->op(2));
 504         bool first_index_tried = false;
 505         int sign = 1;
 506
 507 again:
 508         if (self_idx->is_symbolic()) {
 509                 for (size_t i=1; i<other->nops(); i++) {
 510                         const idx &other_idx = ex_to<idx>(other->op(i));
 511                         if (is_dummy_pair(*self_idx, other_idx)) {
 512
 513                                 // Contraction found, remove metric tensor and substitute
 514                                 // index in second object (assign *self last because this
 515                                 // invalidates free_idx)
 516                                 *other = other->subs(other_idx == *free_idx);
 517                                 *self = (static_cast<const spinidx *>(self_idx)->is_covariant() ? sign : -sign);
 518                                 return true;
 519                         }
 520                 }
 521         }
 522
 523         if (!first_index_tried) {
 524
 525                 // No contraction with first index found, try second index
 526                 self_idx = &ex_to<idx>(self->op(2));
 527                 free_idx = &ex_to<idx>(self->op(1));
 528                 first_index_tried = true;
 529                 sign = -sign;
 530                 goto again;
 531         }
 532
 533         return false;
 534 }
 535
 536 /** Contraction of epsilon tensor with something else. */
 537 bool tensepsilon::contract_with(exvector::iterator self, exvector::iterator other, exvector & v) const
 538 {
 539         GINAC_ASSERT(is_a<indexed>(*self));
 540         GINAC_ASSERT(is_a<indexed>(*other));
 541         GINAC_ASSERT(is_a<tensepsilon>(self->op(0)));
 542         size_t num = self->nops() - 1;
 543
 544         if (is_exactly_a<tensepsilon>(other->op(0)) && num+1 == other->nops()) {
 545
 546                 // Contraction of two epsilon tensors is a determinant
 547                 bool variance = is_a<varidx>(self->op(1));
 548                 matrix M(num, num);
 549                 for (size_t i=0; i<num; i++) {
 550                         for (size_t j=0; j<num; j++) {
 551                                 if (minkowski)
 552                                         M(i, j) = lorentz_g(self->op(i+1), other->op(j+1), pos_sig);
 553                                 else if (variance)
 554                                         M(i, j) = metric_tensor(self->op(i+1), other->op(j+1));
 555                                 else
 556                                         M(i, j) = delta_tensor(self->op(i+1), other->op(j+1));
 557                         }
 558                 }
 559                 int sign = minkowski ? -1 : 1;
 560                 *self = sign * M.determinant().simplify_indexed();
 561                 *other = _ex1;
 562                 return true;
 563         }
 564
 565         return false;
 566 }
 567
 568 //////////
 569 // global functions
 570 //////////
 571
 572 ex delta_tensor(const ex & i1, const ex & i2)
 573 {
 574         static ex delta = (new tensdelta)->setflag(status_flags::dynallocated);
 575
 576         if (!is_a<idx>(i1) || !is_a<idx>(i2))
 577                 throw(std::invalid_argument("indices of delta tensor must be of type idx"));
 578
 579         return indexed(delta, symmetric2(), i1, i2);
 580 }
 581
 582 ex metric_tensor(const ex & i1, const ex & i2)
 583 {
 584         static ex metric = (new tensmetric)->setflag(status_flags::dynallocated);
 585
 586         if (!is_a<varidx>(i1) || !is_a<varidx>(i2))
 587                 throw(std::invalid_argument("indices of metric tensor must be of type varidx"));
 588
 589         return indexed(metric, symmetric2(), i1, i2);
 590 }
 591
 592 ex lorentz_g(const ex & i1, const ex & i2, bool pos_sig)
 593 {
 594         static ex metric_neg = (new minkmetric(false))->setflag(status_flags::dynallocated);
 595         static ex metric_pos = (new minkmetric(true))->setflag(status_flags::dynallocated);
 596
 597         if (!is_a<varidx>(i1) || !is_a<varidx>(i2))
 598                 throw(std::invalid_argument("indices of metric tensor must be of type varidx"));
 599
 600         return indexed(pos_sig ? metric_pos : metric_neg, symmetric2(), i1, i2);
 601 }
 602
 603 ex spinor_metric(const ex & i1, const ex & i2)
 604 {
 605         static ex metric = (new spinmetric)->setflag(status_flags::dynallocated);
 606
 607         if (!is_a<spinidx>(i1) || !is_a<spinidx>(i2))
 608                 throw(std::invalid_argument("indices of spinor metric must be of type spinidx"));
 609         if (!ex_to<idx>(i1).get_dim().is_equal(2) || !ex_to<idx>(i2).get_dim().is_equal(2))
 610                 throw(std::runtime_error("index dimension for spinor metric must be 2"));
 611
 612         return indexed(metric, antisymmetric2(), i1, i2);
 613 }
 614
 615 ex epsilon_tensor(const ex & i1, const ex & i2)
 616 {
 617         static ex epsilon = (new tensepsilon)->setflag(status_flags::dynallocated);
 618
 619         if (!is_a<idx>(i1) || !is_a<idx>(i2))
 620                 throw(std::invalid_argument("indices of epsilon tensor must be of type idx"));
 621
 622         ex dim = ex_to<idx>(i1).get_dim();
 623         if (!dim.is_equal(ex_to<idx>(i2).get_dim()))
 624                 throw(std::invalid_argument("all indices of epsilon tensor must have the same dimension"));
 625         if (!ex_to<idx>(i1).get_dim().is_equal(_ex2))
 626                 throw(std::runtime_error("index dimension of epsilon tensor must match number of indices"));
 627
 628         if(is_a<wildcard>(i1.op(0))||is_a<wildcard>(i2.op(0)))
 629                 return indexed(epsilon, antisymmetric2(), i1, i2).hold();
 630
 631         return indexed(epsilon, antisymmetric2(), i1, i2);
 632 }
 633
 634 ex epsilon_tensor(const ex & i1, const ex & i2, const ex & i3)
 635 {
 636         static ex epsilon = (new tensepsilon)->setflag(status_flags::dynallocated);
 637
 638         if (!is_a<idx>(i1) || !is_a<idx>(i2) || !is_a<idx>(i3))
 639                 throw(std::invalid_argument("indices of epsilon tensor must be of type idx"));
 640
 641         ex dim = ex_to<idx>(i1).get_dim();
 642         if (!dim.is_equal(ex_to<idx>(i2).get_dim()) || !dim.is_equal(ex_to<idx>(i3).get_dim()))
 643                 throw(std::invalid_argument("all indices of epsilon tensor must have the same dimension"));
 644         if (!ex_to<idx>(i1).get_dim().is_equal(_ex3))
 645                 throw(std::runtime_error("index dimension of epsilon tensor must match number of indices"));
 646
 647         if(is_a<wildcard>(i1.op(0))||is_a<wildcard>(i2.op(0))||is_a<wildcard>(i3.op(0)))
 648                 return indexed(epsilon, antisymmetric3(), i1, i2, i3).hold();
 649
 650         return indexed(epsilon, antisymmetric3(), i1, i2, i3);
 651 }
 652
 653 ex lorentz_eps(const ex & i1, const ex & i2, const ex & i3, const ex & i4, bool pos_sig)
 654 {
 655         static ex epsilon_neg = (new tensepsilon(true, false))->setflag(status_flags::dynallocated);
 656         static ex epsilon_pos = (new tensepsilon(true, true))->setflag(status_flags::dynallocated);
 657
 658         if (!is_a<varidx>(i1) || !is_a<varidx>(i2) || !is_a<varidx>(i3) || !is_a<varidx>(i4))
 659                 throw(std::invalid_argument("indices of Lorentz epsilon tensor must be of type varidx"));
 660
 661         ex dim = ex_to<idx>(i1).get_dim();
 662         if (!dim.is_equal(ex_to<idx>(i2).get_dim()) || !dim.is_equal(ex_to<idx>(i3).get_dim()) || !dim.is_equal(ex_to<idx>(i4).get_dim()))
 663                 throw(std::invalid_argument("all indices of epsilon tensor must have the same dimension"));
 664         if (!ex_to<idx>(i1).get_dim().is_equal(_ex4))
 665                 throw(std::runtime_error("index dimension of epsilon tensor must match number of indices"));
 666
 667         if(is_a<wildcard>(i1.op(0))||is_a<wildcard>(i2.op(0))||is_a<wildcard>(i3.op(0))||is_a<wildcard>(i4.op(0)))
 668                 return indexed(pos_sig ? epsilon_pos : epsilon_neg, antisymmetric4(), i1, i2, i3, i4).hold();
 669
 670         return indexed(pos_sig ? epsilon_pos : epsilon_neg, antisymmetric4(), i1, i2, i3, i4);
 671 }
 672
 673 } // namespace GiNaC