ginac/tensor.cpp

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