- archive file format is a little more compact (stores property name and type
[ginac.git] / ginac / lortensor.cpp
1 /** @file lortensor.cpp
2  *
3  *  Implementation of GiNaCĀ“s lortensor objects.
4  *  No real implementation yet, do be done.      */
5
6 /*
7  *  GiNaC Copyright (C) 1999-2000 Johannes Gutenberg University Mainz, Germany
8  *
9  *  This program is free software; you can redistribute it and/or modify
10  *  it under the terms of the GNU General Public License as published by
11  *  the Free Software Foundation; either version 2 of the License, or
12  *  (at your option) any later version.
13  *
14  *  This program is distributed in the hope that it will be useful,
15  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
16  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
17  *  GNU General Public License for more details.
18  *
19  *  You should have received a copy of the GNU General Public License
20  *  along with this program; if not, write to the Free Software
21  *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
22  */
23
24 #include <string>
25 #include <list>
26 #include <algorithm>
27 #include <iostream>
28 #include <stdexcept>
29 #include <map>
30
31 #include "basic.h"
32 #include "add.h"
33 #include "mul.h"
34 #include "debugmsg.h"
35 #include "flags.h"
36 #include "lst.h"
37 #include "lortensor.h"
38 #include "operators.h"
39 #include "tinfos.h"
40 #include "power.h"
41 #include "symbol.h"
42 #include "utils.h"
43 #include "config.h"
44
45 #ifndef NO_NAMESPACE_GINAC
46 namespace GiNaC {
47 #endif // ndef NO_NAMESPACE_GINAC
48
49 //////////
50 // default constructor, destructor, copy constructor assignment operator and helpers
51 //////////
52
53 // public
54
55 lortensor::lortensor()
56 {
57     debugmsg("lortensor default constructor",LOGLEVEL_CONSTRUCT);
58     serial=next_serial++;
59     name=autoname_prefix()+ToString(serial);
60     tinfo_key=TINFO_lortensor;
61 }
62
63 lortensor::~lortensor()
64 {
65     debugmsg("lortensor destructor",LOGLEVEL_DESTRUCT);
66     destroy(0);
67 }
68
69 lortensor::lortensor(const lortensor & other)
70 {
71     debugmsg("lortensor copy constructor",LOGLEVEL_CONSTRUCT);
72     copy (other);
73 }
74
75 const lortensor & lortensor::operator=(const lortensor & other)
76 {
77     debugmsg("lortensor operator=",LOGLEVEL_ASSIGNMENT);
78     if (this != & other) {
79         destroy(1);
80         copy(other);
81     }
82     return *this;
83 }
84
85 //protected
86
87 void lortensor::copy(const lortensor & other)
88 {
89     indexed::copy(other);
90     type=other.type;
91     name=other.name;
92     serial=other.serial;
93 }
94
95 void lortensor::destroy(bool call_parent)
96 {
97     if (call_parent) {
98         indexed::destroy(call_parent);
99     }
100 }
101
102 //////////
103 // other constructors
104 //////////
105
106 // protected
107
108 lortensor::lortensor(lortensor_types const lt, const string & n) : type(lt), name(n)
109 {
110     debugmsg("lortensor constructor from lortensor_types,string",LOGLEVEL_CONSTRUCT);
111     serial=next_serial++;
112     tinfo_key=TINFO_lortensor;
113 }
114
115 lortensor::lortensor(lortensor_types const lt, const string & n, const ex & mu) : indexed(mu), type(lt), name(n)
116 {
117     debugmsg("lortensor constructor from lortensor_types,string,ex",LOGLEVEL_CONSTRUCT);
118     serial=next_serial++;    
119     GINAC_ASSERT(all_of_type_lorentzidx());
120     tinfo_key=TINFO_lortensor;
121 }
122
123 lortensor::lortensor(lortensor_types const lt, const string & n, const ex & mu, const ex & nu) : indexed(mu,nu), type(lt), name(n)
124 {
125     debugmsg("lortensor constructor from lortensor_types,string,ex,ex",LOGLEVEL_CONSTRUCT);
126     serial=next_serial++;
127     GINAC_ASSERT(all_of_type_lorentzidx());
128     tinfo_key=TINFO_lortensor;
129 }
130
131 lortensor::lortensor(lortensor_types const lt, const string & n, const ex & mu, const ex & nu, const ex & rho) : indexed(mu,nu,rho), type(lt), name(n)
132 {
133     debugmsg("lortensor constructor from lortensor_types,string,ex,ex,ex",LOGLEVEL_CONSTRUCT);
134     serial=next_serial++;
135     GINAC_ASSERT(all_of_type_lorentzidx());
136     tinfo_key=TINFO_lortensor;
137 }
138
139 lortensor::lortensor(lortensor_types const lt, const string & n, const ex & mu, const ex & nu, const ex & rho, const ex & sigma) : indexed(mu,nu,rho,sigma), type(lt), name(n)
140 {
141     debugmsg("lortensor constructor from lortensor_types,string,ex,ex,ex,ex",LOGLEVEL_CONSTRUCT);
142     serial=next_serial++;
143     GINAC_ASSERT(all_of_type_lorentzidx());
144     tinfo_key=TINFO_lortensor;
145 }
146
147 lortensor::lortensor(lortensor_types const lt, const string & n, const exvector & iv) : indexed(iv), type(lt), name(n)
148 {
149     debugmsg("lortensor constructor from lortensor_types,string,exvector",LOGLEVEL_CONSTRUCT);
150     serial=next_serial++;
151     GINAC_ASSERT(all_of_type_lorentzidx());
152     tinfo_key=TINFO_lortensor;
153 }
154
155 lortensor::lortensor(lortensor_types const lt, const string & n, unsigned s, const exvector & iv) : indexed(iv), type(lt), name(n), serial(s)
156 {
157     debugmsg("lortensor constructor from lortensor_types,string,unsigned,exvector",LOGLEVEL_CONSTRUCT);
158     GINAC_ASSERT(all_of_type_lorentzidx());
159     tinfo_key=TINFO_lortensor;
160 }
161
162 lortensor::lortensor(lortensor_types const lt, const string & n, unsigned s, exvector *ivp) : indexed(ivp), type(lt), name(n), serial(s)
163 {
164     debugmsg("lortensor constructor from lortensor_types,string,unsigned,exvector",LOGLEVEL_CONSTRUCT);
165     GINAC_ASSERT(all_of_type_lorentzidx());
166     tinfo_key=TINFO_lortensor;
167 }
168
169 //////////
170 // functions overriding virtual functions from bases classes
171 //////////
172
173 //public
174
175 basic * lortensor::duplicate() const
176 {
177     debugmsg("lortensor duplicate",LOGLEVEL_DUPLICATE);
178     return new lortensor(*this);
179 }
180
181 void lortensor::printraw(ostream & os) const
182 {
183     debugmsg("lortensor printraw",LOGLEVEL_PRINT);
184     os << "lortensor(type=" << (unsigned)type
185        << ",indices=";
186     printrawindices(os);
187     os << ",hash=" << hashvalue << ",flags=" << flags << ")";
188 }
189
190 void lortensor::printtree(ostream & os, unsigned indent) const
191 {
192     debugmsg("lortensor printtree",LOGLEVEL_PRINT);
193     os << string(indent,' ') <<"lortensor object: "
194        << "type=" << (unsigned)type << ","
195        << seq.size() << " indices" << endl;
196     printtreeindices(os,indent);
197     os << string(indent,' ') << "hash=" << hashvalue
198        << " (0x" << hex << hashvalue << dec << ")"
199        << ", flags=" << flags << endl;
200 }
201
202 void lortensor::print(ostream & os, unsigned upper_precedence) const
203 {
204     debugmsg("lortensor print",LOGLEVEL_PRINT);
205     switch (type) {
206     case lortensor_g:
207         os << "g";
208         break;
209     case lortensor_rankn:
210         os << name;
211         break;
212     case lortensor_rank1:
213         os << name;
214         break;
215     case lortensor_rank2:
216         os << name;
217         break;
218     case lortensor_epsilon:
219         os << "epsilon";
220         break;
221     case invalid:
222     default:
223         os << "INVALID_LORTENSOR_OBJECT";
224         break;
225     }
226     printindices(os);
227 }
228
229 void lortensor::printcsrc(ostream & os, unsigned type, unsigned upper_precedence) const
230 {
231     debugmsg("lortensor print csrc",LOGLEVEL_PRINT);
232     print(os,upper_precedence);
233 }
234
235 bool lortensor::info(unsigned inf) const
236 {
237     return indexed::info(inf);
238 }
239
240 ex lortensor::eval(int level) const
241 {
242     if (type==lortensor_g) {
243         // canonicalize indices
244         exvector iv=seq;
245         int sig=canonicalize_indices(iv,false); //symmetric
246         if (sig!=INT_MAX) {
247             //something has changed while sorting indices, more evaluations later
248             return ex(sig) *lortensor(type,name,iv);
249         }
250         const lorentzidx & idx1=ex_to_lorentzidx(seq[0]);
251         const lorentzidx & idx2=ex_to_lorentzidx(seq[1]);
252         if ((!idx1.is_symbolic()) && (!idx2.is_symbolic())) {
253             //both indices are numeric
254             if ((idx1.get_value()==idx2.get_value())) {
255                 //both on diagonal
256                 if (idx1.get_value()==0){
257                     // (0,0)
258                     return _ex1();
259                 } else {
260                     if (idx1.is_covariant() != idx2.is_covariant()) {
261                         // (_i,~i) or (~i,_i), i = 1...3
262                         return _ex1();
263                     } else {
264                         // (_i,_i) or (~i,~i), i= 1...3
265                         return _ex_1();
266                     }
267                 }
268             } else {
269                 // at least one off-diagonal
270                 return _ex0();
271             }
272         } else if (idx1.is_symbolic() && idx1.is_co_contra_pair(idx2)) {
273             return Dim()-idx1.get_dim_parallel_space();
274         }
275     }
276     return this -> hold();
277 }
278
279 //protected
280
281 int lortensor::compare_same_type(const basic & other) const
282 {
283     GINAC_ASSERT(is_of_type(other,lortensor));
284     const lortensor *o = static_cast <const lortensor *> (&other);
285     if (type==o->type) {
286         if (type==lortensor_rankn) {
287             if (serial!=o->serial) {
288                 return serial < o->serial ? -1 : 1;
289             }
290         }
291         return indexed::compare_same_type(other);
292     }
293     return type < o->type ? -1 : 1;            
294 }
295
296 bool lortensor::is_equal_same_type(const basic & other) const
297 {
298     GINAC_ASSERT(is_of_type(other,lortensor));
299     const lortensor *o=static_cast<const lortensor *> (&other);
300     if (type!=o->type) return false;
301     if (type==lortensor_rankn) {
302         if (serial!=o->serial) return false;
303     }
304     return indexed::is_equal_same_type(other);            
305 }
306
307 unsigned lortensor::return_type(void) const
308 {
309     return return_types::commutative;
310 }
311 unsigned lortensor::return_type_tinfo(void) const
312 {
313     return tinfo_key;
314 }
315 ex lortensor::thisexprseq(const exvector & v) const
316 {
317     return lortensor(type,name,serial,v);
318 }
319 ex lortensor::thisexprseq(exvector *vp) const
320 {
321     return lortensor(type,name,serial,vp);
322 }
323     
324 //////////
325 // non-virtual functions in this class
326 //////////
327
328 // protected
329
330 void lortensor::setname(const string & n)
331 {
332     name=n;
333 }
334
335 bool lortensor::all_of_type_lorentzidx(void) const
336 {
337     for (exvector::const_iterator cit=seq.begin(); cit!=seq.end(); ++ cit) {
338         if (!is_ex_of_type(*cit,lorentzidx)) {
339             return false;
340         }
341     }
342     return true;
343 }
344
345 // private
346
347 string & lortensor::autoname_prefix(void)
348 {
349     static string * s=new string("lortensor");
350     return *s;
351 }
352
353 //////////
354 // static member variables
355 //////////
356
357 // private
358
359 unsigned lortensor::next_serial=0;
360
361 //////////
362 // friend functions
363 //////////
364
365 lortensor lortensor_g(const ex & mu, const ex & nu)
366 {
367     return lortensor(lortensor::lortensor_g,"",mu,nu);
368 }
369
370 lortensor lortensor_epsilon(const ex & mu, const ex & nu, const ex & rho, const ex & sigma)
371 {
372     return lortensor(lortensor::lortensor_epsilon,"",mu,nu,rho,sigma);
373 }
374
375 lortensor lortensor_rank1(const string & n, const ex & mu)
376 {
377     return lortensor(lortensor::lortensor_rank1,n,mu);
378 }
379
380 lortensor lortensor_rank2(const string & n, const ex & mu, const ex & nu)
381 {
382     return lortensor(lortensor::lortensor_rank2,n,mu,nu);
383 }
384
385 ex simplify_lortensor_mul(const ex & m)
386 {
387     GINAC_ASSERT(is_ex_exactly_of_type(m,mul));
388     exvector v_contracted;
389
390     // collect factors in an exvector, store squares twice
391     int n=m.nops();
392     v_contracted.reserve(2*n);
393     for (int i=0; i<n; ++i) {
394         ex f=m.op(i);
395         if (is_ex_exactly_of_type(f,power)&&f.op(1).is_equal(_ex2())) {
396             v_contracted.push_back(f.op(0));
397             v_contracted.push_back(f.op(0));
398         } else {
399             v_contracted.push_back(f);
400         }
401     }
402
403     unsigned replacements;
404     bool something_changed=false;
405
406     exvector::iterator it=v_contracted.begin();
407     while (it!=v_contracted.end()) {
408         // process only lor_g objects
409         if (is_ex_exactly_of_type(*it,lortensor) &&
410             (ex_to_lortensor(*it).type==lortensor::lortensor_g)) {            
411             const lortensor & g=ex_to_lortensor(*it);
412             GINAC_ASSERT(g.seq.size()==2);
413             const idx & first_idx=ex_to_lorentzidx(g.seq[0]);
414             const idx & second_idx=ex_to_lorentzidx(g.seq[1]);
415             // g_{mu,mu} should have been contracted in lortensor::eval()
416             GINAC_ASSERT(!first_idx.is_equal(second_idx));
417             ex saved_g=*it; // save to restore it later
418
419             // try to contract first index
420             replacements=0;
421             if (first_idx.is_symbolic()) {
422                 replacements = subs_index_in_exvector(v_contracted,
423                                    first_idx.toggle_covariant(),second_idx);
424                 if (replacements==0) {
425                     // not contracted, restore g object
426                     *it=saved_g;
427                 } else {
428                     // a contracted index should occur exactly once
429                     GINAC_ASSERT(replacements==1);
430                     *it=_ex1();
431                     something_changed=true;
432                 }
433             }
434
435             // try second index only if first was not contracted
436             if ((replacements==0)&&(second_idx.is_symbolic())) {
437                 // first index not contracted, *it is again the original g object
438                 replacements = subs_index_in_exvector(v_contracted,
439                                    second_idx.toggle_covariant(),first_idx);
440                 if (replacements==0) {
441                     // not contracted except in itself, restore g object
442                     *it=saved_g;
443                 } else {
444                     // a contracted index should occur exactly once
445                     GINAC_ASSERT(replacements==1);
446                     *it=_ex1();
447                     something_changed=true;
448                 }
449             }
450         }
451         ++it;
452     }
453     if (something_changed) {
454         return mul(v_contracted);
455     }
456     return m;
457 }
458
459 ex simplify_lortensor(const ex & e)
460 {
461     // all simplification is done on expanded objects
462     ex e_expanded=e.expand();
463
464     // simplification of sum=sum of simplifications
465     if (is_ex_exactly_of_type(e_expanded,add)) {
466         ex sum=_ex0();
467         for (unsigned i=0; i<e_expanded.nops(); ++i) {
468             sum += simplify_lortensor(e_expanded.op(i));
469         }
470         return sum;
471     }
472
473     // simplification of commutative product=commutative product of simplifications
474     if (is_ex_exactly_of_type(e_expanded,mul)) {
475         return simplify_lortensor_mul(e);
476     }
477
478     // cannot do anything
479     return e_expanded;
480 }
481
482 ex Dim(void)
483 {
484     static symbol * d=new symbol("dim");
485     return *d;
486 }
487
488 //////////
489 // global constants
490 //////////
491
492 const lortensor some_lortensor;
493 const type_info & typeid_lortensor=typeid(some_lortensor);
494
495 #ifndef NO_NAMESPACE_GINAC
496 } // namespace GiNaC
497 #endif // ndef NO_NAMESPACE_GINAC