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