- g~mu_mu contraction returns Dim for general Lorentz indices and Dim-dimP
[ginac.git] / ginac / lortensor.cpp
1 /** @file lortensor.cpp
2  *
3  *  Implementation of GiNaC's Lorentz 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 <string>
24 #include <list>
25 #include <algorithm>
26 #include <iostream>
27 #include <stdexcept>
28 #include <map>
29
30 #include "basic.h"
31 #include "add.h"
32 #include "mul.h"
33 #include "debugmsg.h"
34 #include "lst.h"
35 #include "lortensor.h"
36 #include "operators.h"
37 #include "tinfos.h"
38 #include "power.h"
39 #include "archive.h"
40 #include "utils.h"
41 #include "config.h"
42
43 #ifndef NO_NAMESPACE_GINAC
44 namespace GiNaC {
45 #endif // ndef NO_NAMESPACE_GINAC
46
47 GINAC_IMPLEMENT_REGISTERED_CLASS(lortensor, indexed)
48
49 //////////
50 // default constructor, destructor, copy constructor assignment operator and helpers
51 //////////
52
53 // public
54
55 lortensor::lortensor() : inherited(TINFO_lortensor), type(invalid)
56 {
57         debugmsg("lortensor default constructor",LOGLEVEL_CONSTRUCT);
58         serial=next_serial++;
59         name=autoname_prefix()+ToString(serial);
60 }
61
62 //protected
63
64 void lortensor::copy(const lortensor & other)
65 {
66         inherited::copy(other);
67         type=other.type;
68         name=other.name;
69         serial=other.serial;
70 }
71
72 void lortensor::destroy(bool call_parent)
73 {
74         if (call_parent) inherited::destroy(call_parent);
75 }
76
77 //////////
78 // other constructors
79 //////////
80
81 // protected
82
83 /** Construct object without any Lorentz index. This constructor is for
84  *  internal use only. */
85 lortensor::lortensor(lortensor_types const lt, const std::string & n) : type(lt), name(n)
86 {
87         debugmsg("lortensor constructor from lortensor_types,string",LOGLEVEL_CONSTRUCT);
88         if (lt == lortensor_symbolic)
89                 serial = next_serial++;
90         else
91                 serial = 0;
92         tinfo_key = TINFO_lortensor;
93 }
94
95 /** Construct object with one Lorentz index. This constructor is for
96  *  internal use only. Use the lortensor_vector() or lortensor_symbolic()
97  *  functions instead.
98  *  @see lortensor_vector
99  *  @see lortensor_symbolic */
100 lortensor::lortensor(lortensor_types const lt, const std::string & n, const ex & mu) : inherited(mu), type(lt), name(n)
101 {
102         debugmsg("lortensor constructor from lortensor_types,string,ex",LOGLEVEL_CONSTRUCT);
103         GINAC_ASSERT(all_of_type_lorentzidx());
104         if (lt == lortensor_symbolic)
105                 serial = next_serial++;
106         else
107                 serial = 0;
108         tinfo_key=TINFO_lortensor;
109 }
110
111 /** Construct object with two Lorentz indices. This constructor is for
112  *  internal use only. Use the lortensor_g(), lortensor_delta() or
113  *  lortensor_symbolic() functions instead.
114  *  @see lortensor_g
115  *  @see lortensor_delta
116  *  @see lortensor_symbolic */
117 lortensor::lortensor(lortensor_types const lt, const std::string & n, const ex & mu, const ex & nu) : inherited(mu,nu), type(lt), name(n)
118 {
119         debugmsg("lortensor constructor from lortensor_types,string,ex,ex",LOGLEVEL_CONSTRUCT);
120         GINAC_ASSERT(all_of_type_lorentzidx());
121         if (lt == lortensor_symbolic)
122                 serial = next_serial++;
123         else
124                 serial = 0;
125         tinfo_key=TINFO_lortensor;
126 }
127
128 /** Construct object with three Lorentz indices. This constructor is for
129  *  internal use only. Use the lortensor_symbolic() function instead.
130  *  @see lortensor_symbolic */
131 lortensor::lortensor(lortensor_types const lt, const std::string & n, const ex & mu, const ex & nu, const ex & rho) : inherited(mu,nu,rho), type(lt), name(n)
132 {
133         debugmsg("lortensor constructor from lortensor_types,string,ex,ex,ex",LOGLEVEL_CONSTRUCT);
134         GINAC_ASSERT(all_of_type_lorentzidx());
135         if (lt == lortensor_symbolic)
136                 serial = next_serial++;
137         else
138                 serial = 0;
139         tinfo_key=TINFO_lortensor;
140 }
141
142 /** Construct object with four Lorentz indices. This constructor is for
143  *  internal use only. Use the lortensor_epsilon() or lortensor_symbolic()
144  *  functions instead.
145  *  @see lortensor_epsilon
146  *  @see lortensor_symbolic */
147 lortensor::lortensor(lortensor_types const lt, const std::string & n, const ex & mu, const ex & nu, const ex & rho, const ex & sigma) : inherited(mu,nu,rho,sigma), type(lt), name(n)
148 {
149         debugmsg("lortensor constructor from lortensor_types,string,ex,ex,ex,ex",LOGLEVEL_CONSTRUCT);
150         GINAC_ASSERT(all_of_type_lorentzidx());
151         if (lt == lortensor_symbolic)
152                 serial = next_serial++;
153         else
154                 serial = 0;
155         tinfo_key=TINFO_lortensor;
156 }
157
158 /** Construct object with arbitrary number of Lorentz indices. This
159  *  constructor is for internal use only. Use the lortensor_symbolic()
160  *  function instead.
161  *
162  *  @see lortensor_symbolic */
163 lortensor::lortensor(lortensor_types const lt, const std::string & n, const exvector & iv) : inherited(iv), type(lt), name(n)
164 {
165         debugmsg("lortensor constructor from lortensor_types,string,exvector",LOGLEVEL_CONSTRUCT);
166         GINAC_ASSERT(all_of_type_lorentzidx());
167         if (lt == lortensor_symbolic)
168                 serial = next_serial++;
169         else
170                 serial = 0;
171         tinfo_key=TINFO_lortensor;
172 }
173
174 lortensor::lortensor(lortensor_types const lt, const std::string & n, unsigned s, const exvector & iv) : indexed(iv), type(lt), name(n), serial(s)
175 {
176         debugmsg("lortensor constructor from lortensor_types,string,unsigned,exvector",LOGLEVEL_CONSTRUCT);
177         GINAC_ASSERT(all_of_type_lorentzidx());
178         tinfo_key=TINFO_lortensor;
179 }
180
181 lortensor::lortensor(lortensor_types const lt, const std::string & n, unsigned s, exvector *ivp) : indexed(ivp), type(lt), name(n), serial(s)
182 {
183         debugmsg("lortensor constructor from lortensor_types,string,unsigned,exvector",LOGLEVEL_CONSTRUCT);
184         GINAC_ASSERT(all_of_type_lorentzidx());
185         tinfo_key=TINFO_lortensor;
186 }
187
188
189 //////////
190 // archiving
191 //////////
192
193 /** Construct object from archive_node. */
194 lortensor::lortensor(const archive_node &n, const lst &sym_lst) : inherited(n, sym_lst)
195 {
196         debugmsg("lortensor constructor from archive_node", LOGLEVEL_CONSTRUCT);
197         unsigned int ty;
198         if (!(n.find_unsigned("type", ty)))
199                 throw (std::runtime_error("unknown lortensor type in archive"));
200         type = (lortensor_types)ty;
201         if (type == lortensor_symbolic) {
202                 serial = next_serial++;
203                 if (!(n.find_string("name", name)))
204                         name = autoname_prefix() + ToString(serial);
205         } else
206                 serial = 0;
207 }
208
209 /** Unarchive the object. */
210 ex lortensor::unarchive(const archive_node &n, const lst &sym_lst)
211 {
212         ex s = (new lortensor(n, sym_lst))->setflag(status_flags::dynallocated);
213
214         if (ex_to_lortensor(s).type == lortensor_symbolic) {
215                 // If lortensor is in sym_lst, return the existing lortensor
216                 for (unsigned i=0; i<sym_lst.nops(); i++) {
217                         if (is_ex_of_type(sym_lst.op(i), lortensor) && (ex_to_lortensor(sym_lst.op(i)).name == ex_to_lortensor(s).name))
218                                 return sym_lst.op(i);
219                 }
220         }
221         return s;
222 }
223
224 /** Archive the object. */
225 void lortensor::archive(archive_node &n) const
226 {
227         inherited::archive(n);
228         n.add_unsigned("type", type);
229         if (type == lortensor_symbolic)
230                 n.add_string("name", name);
231 }
232
233
234 //////////
235 // functions overriding virtual functions from bases classes
236 //////////
237
238 //public
239
240 void lortensor::printraw(std::ostream & os) const
241 {
242         debugmsg("lortensor printraw",LOGLEVEL_PRINT);
243         os << "lortensor(type=" << (unsigned)type
244            << ",indices=";
245         printrawindices(os);
246         os << ",serial=" << serial;
247         os << ",hash=" << hashvalue << ",flags=" << flags << ")";
248 }
249
250 void lortensor::printtree(std::ostream & os, unsigned indent) const
251 {
252         debugmsg("lortensor printtree",LOGLEVEL_PRINT);
253         os << std::string(indent,' ') <<"lortensor object: "
254            << "type=" << (unsigned)type << ","
255            << seq.size() << " indices" << std::endl;
256         printtreeindices(os,indent);
257         os << std::string(indent,' ') << "hash=" << hashvalue
258            << " (0x" << std::hex << hashvalue << std::dec << ")"
259            << ", flags=" << flags << std::endl;
260 }
261
262 void lortensor::print(std::ostream & os, unsigned upper_precedence) const
263 {
264         debugmsg("lortensor print",LOGLEVEL_PRINT);
265         switch (type) {
266         case lortensor_g:
267                 os << "g";
268                 break;
269         case lortensor_delta:
270                 os << "delta";
271                 break;
272         case lortensor_epsilon:
273                 os << "epsilon";
274                 break;
275         case lortensor_symbolic:
276                 os << name;
277                 break;
278         case invalid:
279         default:
280                 os << "INVALID_LORTENSOR_OBJECT";
281                 break;
282         }
283         printindices(os);
284 }
285
286 bool lortensor::info(unsigned inf) const
287 {
288         return inherited::info(inf);
289 }
290
291 ex lortensor::eval(int level) const
292 {
293         if (type==lortensor_g) {
294                 // canonicalize indices
295                 exvector iv=seq;
296                 int sig=canonicalize_indices(iv,false); //symmetric
297                 if (sig!=INT_MAX) {
298                         //something has changed while sorting indices, more evaluations later
299                         return ex(sig) *lortensor(type,name,iv);
300                 }
301                 const lorentzidx & idx1=ex_to_lorentzidx(seq[0]);
302                 const lorentzidx & idx2=ex_to_lorentzidx(seq[1]);
303                 if ((!idx1.is_symbolic()) && (!idx2.is_symbolic())) {
304                         //both indices are numeric
305                         if ((idx1.get_value()==idx2.get_value())) {
306                                 //both on diagonal
307                                 if (idx1.get_value()==0){
308                                         // (0,0)
309                                         return _ex1();
310                                 } else {
311                                         if (idx1.is_covariant() != idx2.is_covariant()) {
312                                                 // (_i,~i) or (~i,_i), i = 1...3
313                                                 return _ex1();
314                                         } else {
315                                                 // (_i,_i) or (~i,~i), i= 1...3
316                                                 return _ex_1();
317                                         }
318                                 }
319                         } else {
320                                 // at least one off-diagonal
321                                 return _ex0();
322                         }
323                 } else if (idx1.is_symbolic() && idx1.is_co_contra_pair(idx2)) {
324                         if (idx1.is_orthogonal_only())
325                                 return Dim() - idx1.get_dim_parallel_space();
326                         else
327                                 return Dim();
328                 }
329         }
330         return this -> hold();
331 }
332
333 //protected
334
335 int lortensor::compare_same_type(const basic & other) const
336 {
337         GINAC_ASSERT(is_of_type(other,lortensor));
338         const lortensor &o = static_cast<const lortensor &>(other);
339
340         if (type!=o.type) {
341                 // different type
342                 return type < o.type ? -1 : 1;
343         }
344
345         if (type == lortensor_symbolic) {
346                 // symbolic, compare serials
347                 if (serial != o.serial) {
348                         return serial < o.serial ? -1 : 1;
349                 }
350         }
351
352         return inherited::compare_same_type(other);
353 }
354
355 bool lortensor::is_equal_same_type(const basic & other) const
356 {
357         GINAC_ASSERT(is_of_type(other,lortensor));
358         const lortensor &o = static_cast<const lortensor &>(other);
359
360         if (type != o.type) return false;
361         if (type == lortensor_symbolic && serial != o.serial) return false;
362         return inherited::is_equal_same_type(other);            
363 }
364
365 unsigned lortensor::return_type(void) const
366 {
367         return return_types::commutative;
368 }
369
370 unsigned lortensor::return_type_tinfo(void) const
371 {
372         return tinfo_key;
373 }
374
375 ex lortensor::thisexprseq(const exvector & v) const
376 {
377         return lortensor(type,name,serial,v);
378 }
379
380 ex lortensor::thisexprseq(exvector *vp) const
381 {
382         return lortensor(type,name,serial,vp);
383 }
384         
385 //////////
386 // non-virtual functions in this class
387 //////////
388
389 // protected
390
391 /** Check whether all indices are of class lorentzidx or a subclass. This
392  *  function is used internally to make sure that all constructed Lorentz
393  *  tensors really carry Lorentz indices and not some other classes. */
394 bool lortensor::all_of_type_lorentzidx(void) const
395 {
396         for (exvector::const_iterator cit=seq.begin(); cit!=seq.end(); ++ cit) {
397                 if (!is_ex_of_type(*cit,lorentzidx)) return false;
398         }
399         return true;
400 }
401
402 // private
403
404 std::string & lortensor::autoname_prefix(void)
405 {
406         static std::string * s = new std::string("lortensor");
407         return *s;
408 }
409
410 //////////
411 // static member variables
412 //////////
413
414 // private
415
416 unsigned lortensor::next_serial=0;
417
418 //////////
419 // friend functions
420 //////////
421
422 /** Construct an object representing the metric tensor g. The indices must
423  *  be of class lorentzidx.
424  *
425  *  @param mu First index
426  *  @param nu Second index
427  *  @return newly constructed object */
428 lortensor lortensor_g(const ex & mu, const ex & nu)
429 {
430         return lortensor(lortensor::lortensor_g,"",mu,nu);
431 }
432
433 /** Construct an object representing the unity matrix delta. The indices
434  *  must be of class lorentzidx.
435  *
436  *  @param mu First index
437  *  @param nu Second index
438  *  @return newly constructed object */
439 lortensor lortensor_delta(const ex & mu, const ex & nu)
440 {
441         return lortensor(lortensor::lortensor_delta,"",mu,nu);
442 }
443
444 /** Construct an object representing the four-dimensional totally
445  *  antisymmetric tensor epsilon. The indices must be of class lorentzidx.
446  *
447  *  @param mu First index
448  *  @param nu Second index
449  *  @param rho Third index
450  *  @param sigma Fourth index
451  *  @return newly constructed object */
452 lortensor lortensor_epsilon(const ex & mu, const ex & nu, const ex & rho, const ex & sigma)
453 {
454         return lortensor(lortensor::lortensor_epsilon,"",mu,nu,rho,sigma);
455 }
456
457 /** Construct an object representing a symbolic Lorentz vector. The index
458  *  must be of class lorentzidx.
459  *
460  *  @param n Symbolic name
461  *  @param mu Index
462  *  @return newly constructed object */
463 lortensor lortensor_vector(const std::string & n, const ex & mu)
464 {
465         return lortensor(lortensor::lortensor_symbolic,n,mu);
466 }
467
468 /** Construct an object representing a symbolic Lorentz tensor of arbitrary
469  *  rank. The indices must be of class lorentzidx.
470  *
471  *  @param n Symbolic name
472  *  @param iv Vector of indices
473  *  @return newly constructed object */
474 lortensor lortensor_symbolic(const std::string & n, const exvector & iv)
475 {
476         return lortensor(lortensor::lortensor_symbolic,n,iv);
477 }
478
479 ex simplify_lortensor_mul(const ex & m)
480 {
481         GINAC_ASSERT(is_ex_exactly_of_type(m,mul));
482         exvector v_contracted;
483
484         // collect factors in an exvector, store squares twice
485         int n=m.nops();
486         v_contracted.reserve(2*n);
487         for (int i=0; i<n; ++i) {
488                 ex f=m.op(i);
489                 if (is_ex_exactly_of_type(f,power)&&f.op(1).is_equal(_ex2())) {
490                         v_contracted.push_back(f.op(0));
491                         v_contracted.push_back(f.op(0));
492                 } else {
493                         v_contracted.push_back(f);
494                 }
495         }
496
497         unsigned replacements;
498         bool something_changed=false;
499
500         exvector::iterator it=v_contracted.begin();
501         while (it!=v_contracted.end()) {
502                 // process only lor_g objects
503                 if (is_ex_exactly_of_type(*it,lortensor) &&
504                         (ex_to_lortensor(*it).type==lortensor::lortensor_g)) {            
505                         const lortensor & g=ex_to_lortensor(*it);
506                         GINAC_ASSERT(g.seq.size()==2);
507                         const idx & first_idx=ex_to_lorentzidx(g.seq[0]);
508                         const idx & second_idx=ex_to_lorentzidx(g.seq[1]);
509                         // g_{mu,mu} should have been contracted in lortensor::eval()
510                         GINAC_ASSERT(!first_idx.is_equal(second_idx));
511                         ex saved_g=*it; // save to restore it later
512
513                         // try to contract first index
514                         replacements=0;
515                         if (first_idx.is_symbolic()) {
516                                 replacements = subs_index_in_exvector(v_contracted,
517                                                                       first_idx.toggle_covariant(),second_idx);
518                                 if (replacements==0) {
519                                         // not contracted, restore g object
520                                         *it=saved_g;
521                                 } else {
522                                         // a contracted index should occur exactly once
523                                         GINAC_ASSERT(replacements==1);
524                                         *it=_ex1();
525                                         something_changed=true;
526                                 }
527                         }
528
529                         // try second index only if first was not contracted
530                         if ((replacements==0)&&(second_idx.is_symbolic())) {
531                                 // first index not contracted, *it is again the original g object
532                                 replacements = subs_index_in_exvector(v_contracted,
533                                                                       second_idx.toggle_covariant(),first_idx);
534                                 if (replacements==0) {
535                                         // not contracted except in itself, restore g object
536                                         *it=saved_g;
537                                 } else {
538                                         // a contracted index should occur exactly once
539                                         GINAC_ASSERT(replacements==1);
540                                         *it=_ex1();
541                                         something_changed=true;
542                                 }
543                         }
544                 }
545                 ++it;
546         }
547         if (something_changed) {
548                 return mul(v_contracted);
549         }
550         return m;
551 }
552
553 /** Perform some simplifications on an expression containing Lorentz tensors. */
554 ex simplify_lortensor(const ex & e)
555 {
556         // all simplification is done on expanded objects
557         ex e_expanded=e.expand();
558
559         // simplification of sum=sum of simplifications
560         if (is_ex_exactly_of_type(e_expanded,add)) {
561                 ex sum=_ex0();
562                 for (unsigned i=0; i<e_expanded.nops(); ++i) {
563                         sum += simplify_lortensor(e_expanded.op(i));
564                 }
565                 return sum;
566         }
567
568         // simplification of (commutative) product
569         if (is_ex_exactly_of_type(e_expanded,mul)) {
570                 return simplify_lortensor_mul(e);
571         }
572
573         // cannot do anything
574         return e_expanded;
575 }
576
577 #ifndef NO_NAMESPACE_GINAC
578 } // namespace GiNaC
579 #endif // ndef NO_NAMESPACE_GINAC