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