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