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