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