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