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