- replaced the Derivative() function by a more resonable fderivative class;
[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 "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